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ABSTRACT 

This  report  reviews  work  done  in  gaining  some  familiarity  with  methods  of 
passive  geolocation,  and  a  search  for  rules  of  thumb  that  might  tell  us  how 
to  optimise  the  geolocation  for  a  given  scenario.  We  first  cover  the  main 
approaches  to  collecting  angle  of  arrival  data  and  point  out  typical  accuracies. 
Following  this  is  an  account  of  the  mathematics  used  to  analyse  this  data  to 
produce  an  estimate  of  an  emitter’s  location.  We  then  give  an  overview  of 
some  of  the  literature,  and  finish  by  demonstrating  a  Matlab  programme  that 
runs  several  geolocation  algorithms.  Simple  rules  of  thumb  that  specify  how  to 
fly  a  baseline  and  take  data,  so  as  to  maximise  the  accuracies  of  the  different 
techniques,  do  not  appear  to  be  readily  derivable. 
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Numerical  Calculations  for  Passive  Geolocation  Scenarios 


EXECUTIVE  SUMMARY 

This  report  has  been  written  with  two  broad  aims.  The  first  is  to  present  an  overview  of 
current  geolocation  techniques  when  bearings-only  data  is  given.  We  give  a  short  overview 
of  geolocation  literature,  as  well  as  covering  the  backgrounds  to  several  techniques;  the 
Cartesian  Pseudo-Linear  Estimator  (which  we  extend  to  handle  the  case  of  a  moving 
emitter),  the  Gauss-Newton  method.  Recursive  Least  Squares,  and  two  types  of  Kalman 
Filter.  These  discussions  are  given  with  a  view  to  being  sufficient  to  enable  the  reader  to 
write  code  that  implements  these  routines,  as  well  as  understanding  some  of  the  relevant 
theory. 

The  second  aim  of  this  report  is  to  analyse  representative  geolocation  scenarios,  in  an 
effort  to  quantify  how  well  each  technique  can  be  expected  to  perform.  We  do  this  by  using 
a  Monte  Carlo  approach  of  repeated  runs  using  software  written  in  Matlab  specifically  for 
this  research.  We  discuss  a  search  for  rules  of  thumb  that  could  direct  how  to  optimise  the 
geolocation  geometry  for  a  scenario  involving  one  emitter  and  a  constant  velocity  receiver. 
In  this  and  other  contexts,  the  Cramer-Rao  bound  is  discussed  in  some  detail.  This  bound 
gives  the  theoretical  limit  on  the  best  accuracy  that  a  geolocation  approach  can  attain. 

The  results  of  analysing  representative  geolocation  scenarios  show  that  no  straightfor¬ 
ward  rules  of  thumb  appear  to  be  derivable;  it  seems  that  nontrivial  geolocation  scenarios 
contain  too  many  parameters  to  allow  each  of  the  techniques’  geolocation  accuracies  to 
be  readily  derivable  using  a  simple  rule.  The  situation  only  gets  worse  for  more  com¬ 
plex  scenarios  involving  multiple  moving  receivers  and  a  single,  perhaps  moving,  emitter. 
Multiple  emitters  are  not  handled  by  this  report;  their  presence  calls  for  more  specialised 
techniques. 
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Glossary,  Conventions,  and  Constants 

CEKF  Cartesian  Extended  Kalman  Filter,  a  geolocation  algorithm. 

CEP  Circular  Error  Probable.  The  radius  of  a  circle  drawn  around  the  emitter  that 
encompasses  50%  of  all  estimates  made  of  its  position.  Thus  50%  of  the  time,  the 
emitter  will  lie  somewhere  within  a  circle  drawn  around  the  estimate.  In  other 
words,  we  can  be  50%  confident  that  the  emitter  lies  within  the  CEP  distance  of 
any  particular  estimate  we  make.  Increasing  this  parameter  to  a  higher  certainty, 
say  95%,  can  increase  the  corresponding  CEP  figure  greatly. 

DF  Direction  Finding. 

DSTO  Defence  Science  Technology  Organisation. 

G— N  Gauss-Newton  geolocation  algorithm. 

MPEKF  Modified  Polar  Extended  Kalman  Filter,  a  geolocation  algorithm. 

CPLE  Cartesian  Pseudo-Linear  Estimator  geolocation  algorithm. 

RLS  Recursive  Least  Squares  geolocation  algorithm. 

TDOA  Time  Difference  Of  Arrival  (also  known  as  DTOA,  Differential  Time  Of  Arrival). 
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1  Overview 


Passive  geolocation  in  a  broad  sense  has  two  rather  separate  parts.  The  starting  point  for 
locating  an  emitter  is  the  direction- finding  process  (DF),  that  is  based  upon  hardware  and 
signal  reception.  Analysis  then  follows  on  the  observed  data;  and  in  a  noisy  world,  this 
draws  heavily  on  a  statistical  approach. 

Direction-finding  techniques  make  up  a  set  of  possibly  complementary  tools  that  are 
used  to  produce  bearing  data  suitable  for  statistical  analysis.  A  general  description  of 
these  is  presented  in  Section  [2]  of  this  report.  Throughout  the  remaining  sections  we 
concentrate  more  on  the  analysis  of  the  bearings  that  the  DF  system  provides,  taking  the 
provision  of  (noisy)  bearings  as  a  given.  Analysis  of  this  bearing  data  itself  is  described  in 
Section  ^  There  we  discuss  several  different  methods,  all  of  which  are  aimed  at  producing 
an  estimate  of  where  the  emitter  actually  is  and  how  fast  it  might  be  moving. 

Although  geolocation  has  a  long  history,  much  good  work  in  the  held  has  been  done 
only  relatively  recently,  as  computer  power  continues  to  increase.  The  various  techniques 
that  in  the  past  were  graded  on  their  computational  intensity  are  now  all  performed  easily 
at  high  speed  on  a  modern  computer,  and  the  question  of  which  one  is  to  be  preferred  is 
partly  a  question  of  how  easy  they  each  are  to  implement  and  understand.  We  discuss 
recent  work  in  the  various  geolocation  techniques  in  Section  @1 

Section  [5]  deals  with  a  very  specihc  scenario  where  three  hxed  receivers  geolocate  one 
stationary  emitter.  We  produce  Cramer-Rao  bounds,  as  well  as  comparing  the  results  of 
different  geolocation  routines.  Long-baseline  TDOA  is  then  discussed  for  this  geometry, 
with  results  presented  that  give  an  idea  of  how  the  geolocation  is  affected  by  baseline 
length  and  timing  accuracy. 

Since  one  of  the  aims  of  the  current  work  has  been  to  develop  a  feel  for  the  types  of 
paths  that  must  be  flown  to  achieve  a  given  geolocation  accuracy,  we  need  to  be  able  to 
run  the  various  bearing  analysis  routines  on  sample  data  to  see  what  effects  the  various 
parameters  have  on  the  final  estimate  of  the  emitter’s  position  and  state  of  motion.  To  this 
end,  in  Appendix  [C]  we  show  the  results  of  a  Matlab  programme  that  runs  these  routines 
as  required.  Section  [6]  analyses  these  results  with  a  view  to  obtaining  some  rules  of  thumb 
for  simple  scenarios. 


2  Direction-Finding  Techniques 

The  art  of  direction  finding  itself  goes  back  to  Hertz  and  the  beginnings  of  radio  theory, 
when  it  was  realised  that  the  anisotropy  in  an  antenna’s  reception  pattern  could  be  used 
to  hnd  an  emitter’s  direction.  This  is  still  a  common  means  of  direction  Ending,  but  has 
long  been  supplemented  by  more  sophisticated  coherent  techniques.  Because  of  this,  it  is 
convenient  to  split  DF  apparatus  into  two  classes:  those  that  make  use  of  the  receiver’s 
amplitude  response,  and  those  that  use  its  phase  response.  The  amplitude  response  type 
can  be  of  small  size  and  thus  fit  well  into  fighter  aircraft;  on  the  other  hand,  the  accurate 
measurement  of  phase  differences  in  the  arrival  of  a  wave  front  can  require  a  larger  receiver, 
and  antenna  array,  than  might  be  practical  on  these. 
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Even  if  we  are  able  to  locate  an  emitter  to  any  desired  accuracy,  our  measurements 
can  still  be  limited  by  our  knowledge  of  our  own  ship’s  position  and  heading.  However, 
own-ship  position  and  heading  are  determined  very  accurately  using  GPS/INS. 


2.1  Amplitude  Response 

The  hrst  of  the  noncoherent  DF  techniques  is  the  simplest  and  oldest:  locating  the  di¬ 
rection  of  maximum  signal  strength  while  rotating  the  antenna.  This  has  the  simplifying 
advantage  that  the  only  detailed  knowledge  of  the  beam  pattern  needed  is  the  beam  width, 
where  a  typical  resolution  of  this  type  of  system  is  about  1/10  of  that  width.  But  this 
simplicity  brings  several  possible  problems: 

•  The  receiver  might  have  to  be  held  wide  open  to  maximise  signal  reception  probabil¬ 
ity,  which  is  a  disadvantage  in  the  presence  of  many  signals.  On  the  other  hand,  this 
sort  of  wide  band  DF  is  not  always  required,  with  narrow  band  searching  instead 
being  favoured.  In  some  situations  a  wide  band  scan  might  still  be  done,  the  signal 
located  and  then  further  searching  done  over  a  much  smaller  bandwidth. 

•  Many  pulses  need  to  be  intercepted — with  a  correspondingly  large  amount  of  data 
processed. 

•  Signals  of  varying  strength  can  completely  distort  the  measurement.  This  is  the 
main  problem  and  can  lead  to  the  method  being  effectively  useless,  unless  a  cor¬ 
rection  is  applied.  Normally  this  is  overcome  by  employing  an  extra  nonscanning 
nondirectional  antenna  that  also  processes  the  pulses,  allowing  us  to  normalise  the 
data  by  comparing  the  ratios  of  the  signal  strengths  seen  by  each  antenna. 

Alternatively,  we  can  dispense  with  the  rotation  by  using  at  least  two  stationary  antennas 
with  known,  nonisotropic  beam  patterns:  overlapping  these  patterns  and  comparing  the 
signal  strength  from  each  then  gives  the  source  direction,  even  when  the  signal  strength 
is  changing.  This  method  of  amplitude  comparison  is  used  very  commonly  for  direction 
finding.  Typically  there  are  four  or  more  antennas,  with  the  coarse  tuning  being  done  by 
simply  measuring  which  has  the  strongest  signal,  followed  by  hne  tuning  by  comparing 
responses  as  above.  Small  errors  in  power  received  can  lead  to  a  large  DF  error,  so  that 
typical  bearing  accuracy  for  four  receivers  is  10-15°.  This  can  be  increased  to  about  5° 
for  six  receivers,  while  2°  is  attainable  with  eight  receivers. 


2.2  Time  Difference  of  Arrival 

If  our  platform  is  large  enough,  measuring  the  difference  in  arrival  times  of  the  signal 
at  two  or  more  widely  spaced  receivers  will  give  good  DF  information.  This  technique  is 
known  as  Time  Difference  of  Arrival  (TDOA).  In  two  dimensions,  for  each  pair  of  antennas, 
the  emitter  will  be  located  somewhere  on  a  hyperbola;  so  a  third  antenna  on  a  different 
baseline  will  narrow  the  position  down  to  a  point.  In  three  dimensions  the  hyperbolae 
become  surfaces  of  revolution  about  the  line  joining  each  pair,  so  four  receivers  will  be 
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needed.  These  four,  or  commonly  still  more  receivers,  will  be  set  in  a  number  of  different 
baselines. 

For  long  baselines  (greater  than  a  few  tens  of  metres),  this  is  a  modification  to  the 
meaning  of  direction  finding,  in  the  sense  that  we  are  not  locating  the  emitter  on  a  ray 
described  by  one  bearing.  Rather,  a  hyperbola  needs  to  be  specified.  The  situation  is 
simplified  for  baselines  of  perhaps  ten  or  twenty  metres  (“short-baseline  TDOA”),  since 
the  multiple  bearings  are  then  essentially  the  same,  lying  more  or  less  on  the  hyperbola’s 
asymptotes.  This  produces  a  straight  line  geometry  and  a  simple  direction.  In  this  case 
the  receivers  can  be  sited  on  one  platform,  but  the  timing  accuracy  required  then  becomes 
crucial.  The  hyperbolic  case  becomes  more  important,  and  in  principle  easier  to  imple¬ 
ment,  when  the  receivers  are  more  widely  separated.  But  in  practice,  widely  separated 
receivers  are  limited  by  the  requirement  that  they  both  see  the  same  pulse,  and  this  is  to 
some  extent  dictated  by  the  emitter  beamwidth.  This  topic  is  covered  in  some  detail  in 
Section  15.21  as  well  as  in  [1] .  Typical  bearing  accuracy  is  1-2° . 


2.3  Interferometry 

Whereas  TDOA  measures  the  difference  in  arrival  times  of  a  wave  at  two  receivers,  inter¬ 
ferometry  measures  the  difference  in  arrival  phase  at  two  receivers. 

Beasley  and  Miles  [2]  has  a  good  discussion  on  the  various  sources  of  error  possible. 
Notable  errors  result  from  the  following: 

•  Multipath  interference  due  to  the  ESM  platform’s  fuselage  etc.  This  might  produce 
bearings  with  significant  errors,  but  if  the  receivers  are  changing  orientation  from 
one  measurement  to  the  next  (e.g.  if  the  platform  is  turning),  then  these  errors  will 
probably  not  lead  to  any  bias  in  the  measurements. 

•  Multipath  interference  caused  by  the  radome,  which  is  ever-present,  but  is  at  least 
a  systematic  error. 

•  Bearing  ambiguity.  This  will  not  occur  if  the  antenna  spacing  is  less  than  half  a 
wavelength.  Typical  radar  wavelengths  of  3-30  cm  restrict  the  required  antenna 
spacing,  but  in  practice  ambiguities  are  overcome  by  using  multiple  baselines. 

•  Thermal  noise.  This  cannot  easily  be  removed  through  system  design,  and  hence 
may  well  be  the  greatest  error  contributor. 

Interferometry  tends  to  give  very  high  accuracy;  better  than  0.5°  is  possible  for  single 
pulses  on  line-of-sight  paths. 


2.4  Absolute  Doppler 


When  an  emitter  is  at  rest — for  example  anchored  to  Earth — its  Doppler  shift  will  give  us 
information  about  where  it  is.  This  can  be  seen  by  realising  that  provided  we  know  the 
emitter’s  velocity  relative  to  us,  together  with  its  velocity  radially  from  us,  then  we  can 
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in  essence  shift  that  velocity  vector  around  our  own  position,  until  its  radial  component 
matches  the  measured  value.  All  such  matches  can  only  lie  on  a  cone  extending  out  from 
ourselves,  and  this  cone  cuts  Earth’s  surface  in  a  hyperbola.  So  we’ll  know  that  the  emitter 
lies  somewhere  on  this  hyperbola. 

We  do  know  the  emitter’s  velocity  relative  to  us  because  it  is  anchored  to  Earth,  and 
so  is  just  equal  and  opposite  to  our  known  velocity  relative  to  Earth.  We  also  know  its 
radial  velocity,  since  this  is  given  by  the  ratio  of  observed  to  emitted  frequencies: 

fo  ^radial  /n  i  \ 

s"'-— ■ 

The  only  real  problem  here  is  that  we  need  to  know  the  emitted  frequency.  And  in  practice, 
accurate  knowledge  of  this  is  probably  not  available — especially  because  this  frequency  can 
change  very  quickly  in  modern  radars.  Thus,  the  method  might  only  be  of  use  for  narrow- 
band  signals  produced  by  crystal-controlled  radar,  instead  of  the  more  common  magnetron 
types  whose  frequencies  are  much  less  constant. 


2.5  Frequency  Difference  of  Arrival 

Again  for  a  stationary  emitter.  Frequency  Difference  of  Arrival  (FDOA)  uses  the  difference 
in  Doppler  shifts  in  emitter  frequency  as  observed  by  at  least  two  platforms  that  both  move. 
Since  the  Doppler  shifts  yield  the  radial  velocity  of  each  receiver  from  the  emitter,  it’s  a 
straightforward  matter  to  determine  the  emitter  location,  which  will  be  unique  (or  perhaps 
two-fold  ambiguous)  for  most  receiver  configurations. 


2.6  Differential  Doppler 

The  Differential  Doppler  technique  follows  the  same  basic  idea  as  FDOA.  But  whereas 
FDOA  makes  use  of  the  same  signal  received  by  two  receivers  at  widely  separated  points 
in  space,  in  principle  Differential  Doppler  needs  just  one  receiver  to  measure  a  change  in 
frequency  over  time,  as  the  emitter  moves  relative  to  it.  The  complication  arising  is  that 
the  emitted  frequency  might  change  over  time  in  an  unknown  way. 

Beasley  and  Miles  [2]  discuss  Differential  Doppler,  but  actually  have  not  properly 
incorporated  the  differential  Doppler  shift  into  their  analysis.  What  they  are  calculating, 
the  rate  of  increase  of  interferometer  phase  difference  with  time,  has  two  parts:  a  pure 
Doppler  component  (i.e.  involving  dA/dt)  which  they  don’t  include  in  their  analysis,  and 
a  sort  of  FDOA  component  caused  by  the  motion  of  the  emitter  as  seen  by  the  receiver. 
This  last  term  is  the  only  one  written  down  by  Beasley  and  Miles,  although  it  turns  out 
to  be  the  dominant  one  anyway.  They  doubt  this  technique  has  been  tested.  (A  change  is 
necessary  to  their  equation  (2-22):  replace  sin^  (/>  with  cos^  (j)  since  they  have  confused  the 
two  (p’s  in  their  Figures  2-1  and  2-8.) 
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North 


Figure  1:  Bearings  taken  at  the  and  k+1^^  reeeiver  positions 


of  signal.  With  wide-band  signals  the  situation  is  reversed:  the  spread  of  frequencies 
means  Differential  Doppler  does  not  perform  well,  so  that  TDOA  is  then  preferred. 


3  Statistical  Processing  of  Bearing  Data 

3.1  Layout  of  a  Geolocation  Problem 

For  the  next  subsections,  refer  to  the  layout  of  Figure  H  adapted  from  [3].  A  stationary 
emitter  is  situated  at  {xe,ye),  while  a  receiver  moving  along  some  path  takes  bearing 
measurements  at  various  points.  The  receivers  are  situated  at  ri . .  .r„,  and  take  (noisy) 
bearings  /3i . .  ./?„,  which  form  the  data  set.  What  we  wish  to  do  is  either  produce  an 
estimate  of  the  emitter’s  position  (xe,  Ve),  or  else,  given  such  an  estimate,  improve  upon  it 
by  producing  a  new  estimate.  We  also  require  the  error  in  the  result.  In  general  this  can 
be  difficult  to  calculate,  but  we  can  be  guided  as  to  what  error  can  reasonably  be  expected 
by  considering  the  Cramer-Rao  lower  bound  for  the  given  scenario,  as  well  as  the  circular 
error  probable  (CEP)  and  range  errors,  as  discussed  next. 


3.2  Cramer-Rao  Lower  Bound 

The  Cramer-Rao  lower  bound  is  an  indicator  of  how  well  we  can  locate  a  parameter  that 
is,  in  some  sense,  the  “centre”  of  a  distribution.  Its  use  is  motivated  by  this  fact,  since  it 
will  tell  us  how  well  a  geolocation  algorithm  is  performing  compared  to  the  absolute  best 
possible  performance  theoretically  attainable.  This  absolute  benchmark  also  allows  for  a 
better  comparison  of  the  various  algorithms,  rather  than  just  comparing  them  with  each 
other. 
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Suppose  we  receive  a  signal  from  which  we  wish  to  extract  an  interesting  parameter  x. 
For  example,  in  the  case  of  a  normal  distribution,  x  might  be  the  mean.  In  general,  if  a 
data  set  depends  on  this  unknown  parameter,  then  in  the  absence  of  our  being  able  to 
determine  x,  we  wish  to  define  an  “estimator”  x  of  x  in  such  a  way  that  x  is  unbiased;  by 
this  is  meant  that  calculating  it  from  a  collection  of  data  sets  would  yield  x  on  the  average, 
so  that  the  expected  value  E{x}  equals  x.  In  some  situations  the  Cramer-Rao  theory  will 
produce  this  unbiased  estimator  for  us,  as  well  as  a  lower  bound  on  the  associated  variance. 

In  an  unpublished  paper  by  Kim  Brown  at  DSTO  [4],  this  theory  has  been  applied 
to  the  case  of  an  aircraft  flying  at  various  broadside  angles  and  ranges  to  an  emitter,  in 
an  effort  to  calculate  just  what  the  minimum  geolocation  error  will  be,  for  the  case  of 
a  physically  realistic  probability  density  function  of  bearing  errors  seen  by  the  receiver. 
This  function  is  chosen  in  [4j  to  be  von  Mises,  a  symmetrical  bell-shaped  curve  with  its 
maximum  at  zero  degrees  error,  being  the  angular  version  of  a  gaussian  distribution.  (That 
is,  it  gives  the  probability  density  that  the  bearing  suffers  from  some  error,  which  can  be 
positive  or  negative  corresponding  to  the  measured  direction  being  too  far  off  true  in  one 
direction  or  the  other.  The  width  of  the  bell  curve  is,  as  usual,  related  to  the  actual 
bearing  error  we  attribute  to  the  DF  equipment.) 

The  results  in  [4]  have  been  reproduced  independently  using  Matlab  in  Appendix  O 
The  Matlab  programme  goes  further  than  [4],  in  the  sense  that  it  allows  for  arbitrary 
broadside  angles,  emitter  ranges,  receiver  speeds,  observation  periods  and  number  of  stan¬ 
dard  deviations  for  the  uncertainty  ellipse,  as  opposed  to  the  paper’s  static  values.  The 
programme  also  does  something  a  little  different  to  the  paper:  it  applies  the  von  Mises 
probability  distribution  to  the  geolocation  scenarios  it  is  analysing  (as  opposed  to  the  ex¬ 
ample  in  [4j  which  does  not  incorporate  any  tracking  algorithms),  and  plots  an  uncertainty 
ellipse  superimposed  on  the  map  with  the  tracking  algorithm  output. 

This  uncertainty  ellipse  gives  us  an  idea  of  the  best  accuracy  we  can  hope  for  when 
processing  DF  data,  regardless  of  which  method  is  being  used  to  locate  the  emitter  (al¬ 
ways  assuming  a  von  Mises  error,  which  is  quite  reasonable).  The  fact  that  the  parameters 
can  be  freely  input  to  Matlab  means  that  we  can  use  the  Cramer-Rao  theory  as  a  ba¬ 
sis  for  comparing  the  techniques  of  Cartesian  Pseudo-Linear  Estimator,  Gauss-Newton, 
Recursive  Least  Squares,  and  Kalman  filtering,  as  will  be  discussed  shortly. 

Cramer-Rao  bounds  are  calculated  for  a  multiple  receiver  scenario  in  Section  5.1i  It 
must  be  added  that  the  lower  bound  is,  strictly  speaking,  only  applicable  to  an  unbiased 
estimator.  But  it  turns  out  that  the  results  of  geolocation  algorithms  will  generally  tend 
to  be  biased,  so  that  we  must  treat  Cramer-Rao  results  cautiously. 

3.2.1  Two  Types  of  Cramer-Rao  Bound 

Signal  processing  theory  involves  two  types  of  Cramer-Rao  bound.  These  correspond  to 
whether  we  use  non-bayesian  statistics  or  bayesian  statistics  for  the  geolocation  calculation. 
In  the  sections  that  follow,  we  will  describe  several  geolocation  methods.  The  Cartesian 
Pseudo-Linear  Estimator  and  Gauss-Newton  are  “batch”  methods  where  the  state  of  the 
emitter  (position,  velocity,  etc.)  is  treated  as  having  a  definite  but  unknown  value,  with 
no  prior  knowledge  assumed:  this  is  a  non-bayesian  approach.  One  application  of  the 
batch  method  is  sufficient  to  yield  an  estimate  of  that  state.  The  last  three  methods 
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are  all  “recursive” ,  in  which  the  state  is  treated  as  a  random  variable,  according  to  some 
probability  density  function  given  initially  by  prior  knowledge.  This  prior  knowledge 
indicates  that  the  Cramer-Rao  calculation  should  assume  a  bayesian  approach,  since  Bayes 
theory  makes  use  of  this  prior  knowledge. 

The  lower  bound  on  the  variance  of  the  unbiased  estimator  can  be  shown  to  be  the 
inverse  of  the  Fisher  information  J  (a  matrix,  in  general).  That  is,  if  our  unbiased 
estimator  of  the  emitter’s  state  is  x,  then  Cramer-Rao  theory  states  that  the  following  is 
always  true: 

var(x)  ^  J~^  =  CRLB,  the  Cramer-Rao  lower  bound,  (3.1) 

where  the  meaning  of  “greater  than  a  matrix”  is  quantified  just  after  (3.6)  ahead.  How 
do  we  calculate  the  Fisher  information  matrix?  Consider  first  the  non-bayesian  case  that 
is  applied  to  batch  processes. 


Non-Bayesian  Case.  Suppose  the  data  set  of  k  measurements  is 

Zk  =  {z,,...z^},  (3.2) 

which  depends  on  an  unknown  parameter  x  (e.g.,  z^  are  bearing/range  pairs,  and  x  is  the 
position  of  the  emitter).  The  set  might  also  contain  a  first  element  Zq,  which  is  not  a 
measurement  as  such,  but  rather  is  prior  information.  The  probability  that  given  some 
particular  x,  the  k  measurements  Zk  will  be  observed  is  p{Zk\x),  called  the  likelihood.  R 
turns  out  to  be  more  convenient  to  work  with  the  negative  logarithm  of  the  likelihood: 

L{Zk\x)  =  -  \np{Zk\x) .  (3.3) 

For  a  one-dimensional  case  where  x  is  a  scalar,  the  non-bayesian  Fisher  information 
(a  scalar)  is  defined  as  an  expected  value  over  the  data: 

Jnb  ^  Ez.  =  Ez.  I  (|1)  I  (3,4) 


(where  the  last  equality  can  be  proved  easily) .  Generally  the  data  is  dependent  on  several 
parameters,  whose  values  at  the  latest  time  k  are  collectively  called  x^.  Dropping  the 
k  subscripts,  the  Fisher  information  at  time  k  is  now  a  matrix^ 


Jnb  {Va.  {ViL{Z\x))}  =Ez  {{ViL)  {V^L)}  ,  (3.5) 


where  VL  is  a  row  vector,  and  V*L  =  iVL)^  is  a  column  vector.  Component- wise,  the 
previous  expression  is: 


\d‘^L{z\x)\ 

\  dxi  dxj  j  ^  \  dxi  dxj  j 


(3.6) 


(remembering  that  Xi,  xj  are  really  the  i,  j  components  of  Xk).  The  inverse  Fisher  informa¬ 
tion  is  a  lower  bound  in  the  sense  that  cov{x)  —  is  positive  semidefinite — meaning 
that  the  quadratic  form  it  produces  is  always  ^  0. 

^In  general,  the  CRLB  theory  assumes  the  regularity  condition  to  hold:  Ez^{Vxf,L}  =  0  for  all  Xk, 
which  is  true  for  well  behaved  distributions. 
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Bayesian  Case.  The  bayesian  viewpoint  used  to  calculate  the  Cramer~Rao  bound, 
known  as  the  “posterior”  CRLB,  for  recursive  routines  is  much  more  involved,  and  has  been 
briefly  described  in  Appendix  El  It  uses  some  recursive  notation  laid  out  in  Section  13.101 
Good  references  to  its  discussion  of  the  posterior  CRLB  are  [51  [6] . 


3.3  Circular  Error  Probable  and  Range  Errors  as 
Ways  of  Quantifying  Error 

There  are  various  ways  to  quantify  the  accuracy  of  a  geolocation  method.  Two  in  particular 
are  used  in  this  report.  The  first  is  the  circular  error  probable  (CEP).  This  is  the  radius 
of  a  circle  that,  when  drawn  around  the  actual  position  of  an  emitter,  will  encompass 
half  of  a  large  number  of  estimates  to  that  position.  So  this  circle  comprises  the  region 
within  which  there  is  a  50%  probability  of  estimating  the  emitter  to  lie;  hence  the  smaller 
the  CEP,  the  better  is  the  geolocation  method.  Broadly  speaking,  since  the  area  of  the 
circle  is  proportional  to  the  square  of  the  CEP,  a  halving  of  the  CEP  really  signals  a 
dramatic  improvement  in  accuracy.  Although  the  collection  of  estimates  made  to  the 
emitter  position  will  not,  in  general,  lie  symmetrically  around  the  true  emitter  position, 
nevertheless,  the  CEP  is  a  good  one-parameter  quantifier  of  the  accuracy  of  a  geolocation 
method.  In  practice,  the  CEP  will  be  drawn  around  the  estimated  position  of  the  emitter, 
because  the  actual  emitter  position  is  unknown. 

The  second  way  to  quantify  the  accuracy  of  a  geolocation  method  is  to  specify  two 
numbers:  the  cross-  and  down-range  errors.  Cross-range  error  is  the  transverse  difference 
between  the  actual  emitter  position  and  our  estimate  of  it.  Down-range  error  is  the 
difference  between  the  range  to  the  actual  emitter  and  the  range  of  the  estimate,  which  is 
more  or  less  as  depicted  in  Figure  [2]  if  the  cross-range  error  is  not  too  large.  The  figure 
shows  a  set  of  estimates  to  an  emitter’s  position,  with  the  CEP  and  range  errors  indicated. 


3.4  Observability  of  the  Emitter 

In  order  to  determine  the  receiver  network  necessary  for  a  given  geolocation  problem — or 
perhaps  the  parameters  of  the  flight  path  of  one  receiver — we  need  to  give  some  consid¬ 
eration  to  whether  an  emitter’s  position  can  be  uniquely  determined,  even  in  the  absence 
of  noise.  Since  our  measurements  are  bearings  only,  we  can  imagine  a  scenario  in  which 
another,  more  agile,  emitter  moves  arbitrarily,  but  always  so  as  to  remain  hidden  behind 
the  real  emitter.  Thus  the  position  of  an  emitter  is  never  uniquely  determined,  since  this 
other  more  agile  emitter  could  really  have  been  the  only  one  present. 

But  this  more  agile  emitter  must  have  to  follow  a  much  more  convoluted  course  than 
the  visible  one,  and  the  wisest  choice  might  be  to  say  that  whichever  of  these  alternative 
emitters  is  moving  in  the  simplest  way  is  the  real  one.  For  example,  suppose  the  receiver  is 
moving  at  constant  velocity:  that  is,  a  straight  line  with  constant  speed.  Then  an  emitter 
that  is  really  stationary  will  give  rise  to  the  same  bearings  as  one  that  moves  at  just  the 
right  speed  to  always  appear  to  be  behind  or  in  front  of  the  the  real  one.  And  these  could 
hide  another  one  that  is  zig-zagging  wildly  back  and  forth,  again  always  so  as  to  remain 
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Figure  2:  Definitions  of  CEP  and  range  errors 


hidden  behind  the  real  one.  Our  intuition  tells  us  that  those  two  that  move  are  almost 
certainly  not  real,  but  rather  are  “ghosts” .  Depending  on  what  prior  information  we  have, 
the  most  pragmatic  decision  might  be  that  the  real  emitter  is  in  fact  stationary. 

Accepting  the  fact  that  the  real  emitter  could  be  moving  quite  wildly  in  just  the  right 
way  so  as  to  give  rise  to  the  observed  bearings,  we  denote  the  stationary  emitter  in  the 
above  example  observable  if,  when  we  model  it  as  stationary,  a  geolocation  analysis  will 
produce  an  accurate  value  for  its  position.  Likewise  a  constant-velocity  emitter  is  denoted 
as  observable  if,  when  it’s  modelled  as  having  a  constant  velocity,  a  geolocation  analysis 
returns  an  accurate  value  for  its  position  and  that  constant  velocity. 

Relevant  to  this  discussion  is  the  idea  of  manoeuvring,  which  can  be  quantihed  in  the 
following  way.  If  two  objects  A  and  B  are  moving,  then  A  manoeuvres  more  than  B  if, 
for  example,  A  has  a  constant  velocity  while  B  only  has  a  constant  position;  or  A  has  a 
constant  acceleration  while  B  only  has  at  most  a  constant  velocity,  and  so  on. 

It  can  be  shown  that  if  the  emitter  is  observed,  then  the  receiver  must  have  manoeuvred 
more  to  have  found  it.  The  converse  is  not  quite  true:  just  being  more  manoeuvrable  does 
not  guarantee  that  we  will  locate  the  emitter,  but  we  usually  will  even  so,  provided  we 
have  made  a  reasonable  model  of  its  motion.  The  bottom  line  is  that  if  we  (the  receiver) 
suspect  the  emitter  is  stationary,  then  we  must  at  least  be  moving  to  geolocate  it;  while  if 
we  suspect  the  emitter  has  constant  velocity  then  we  must  accelerate,  even  if  that  is  limited 
to  a  single  turn  connecting  two  constant-velocity  paths.  While  an  acceptable  acceleration 
might  be  nothing  more  than  a  kink  in  our  otherwise  constant-velocity  flight  path,  the  fact 
that  this  kink  might  be  very  slight  shows  that  there  is  a  continuum  of  degrees  to  which  an 
emitter  is  observable.  It  is  not  true  that  emitters  are  either  observable  or  not;  the  accuracy 
of  a  geolocation  analysis  depends  on  how  we  model  the  emitter’s  motion,  and  how  we  (the 
receiver)  move.  Different  geolocation  techniques  have  differing  degrees  of  sensitivity  to  the 
emitter’s  observability. 
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All  of  this  has  the  consequence  of  limiting  what  can  be  achieved,  even  in  principle,  with 
target  motion  analysis.  When  implementing  a  geolocation  routine,  we  must  specify  the 
type  of  emitter  motion:  for  example,  if  we  suspect  that  the  emitter  has  constant  velocity, 
then  we  (the  receiver)  must  accelerate  to  out-manoeuvre  it  in  order  for  the  geolocation  to 
have  any  chance  of  making  sense;  and  we  also  will  need  to  provide  some  sort  of  model  of 
the  emitter  that  allows  it  to  have  a  constant  velocity. 

It  might  be  thought  that  since  constant  velocity  is  only  a  special  kind  of  constant 
acceleration  (where  the  acceleration  is  just  zero),  it  might  be  useful  to  model  the  emitter 
as  accelerating,  even  though  we  suspect  it  probably  only  has  a  constant  velocity.  Or,  to 
take  the  situation  to  its  logical  conclusion,  we  might  be  able  to  model  the  emitter  as  having 
some  high  manoeuvrability,  and  so  long  as  the  receiver  outmanoeuvres  it  in  actuality,  then 
it  should  be  able  to  find  the  emitter’s  true  motion.  But  this  cannot  be  true.  Since  we  can 
imagine  that  there  is  an  infinite  hierarchy  of  ever  more  complicated  manoeuvring  emitters 
hidden  behind  the  real  one,  as  soon  as  we  model  the  emitter  as  anything  of  a  high  order 
motion,  we  will  open  ourselves  up  to  confusing  its  motion  with  all  of  the  ghost  emitters 
that  have  manoeuvrability  at  most  as  complicated  as  that  of  our  model.  So  there  will  be 
an  ambiguity  reflected  in  our  results  of  where  the  emitter  is. 

In  fact  in  practice,  for  some  geolocation  techniques  (“smoothers”  such  as  the  CPLE 
technique  described  in  Sect.  3.7).  the  more  manoeuvrability  we  build  into  the  model  of 
the  emitter,  the  closer  its  estimated  position  and  motion  lie  to  the  receiver’s  position  and 
motion.  It  actually  appears  to  be  sitting  right  on  top  of  the  receiver.  With  hindsight  this 
is  not  as  strange  a  result  as  it  sounds,  because  such  an  emitter  will  certainly  always  be 
sighted  in  the  direction  given  by  the  measured  bearings.  It’s  as  if  the  geolocation  routine 
is  giving  up,  and  simply  returning  an  emitter  state  that  will  always  “fit”,  in  the  sense  of 
being  compatible  with  the  measured  bearings. 


3.5  Producing  an  Estimate  of  the  Emitter  Position: 

The  Stansfield  Algorithm 

Perhaps  the  earliest  of  the  “modern”  papers  to  discuss  the  processing  of  bearing  data  was 
published  by  Stansfield  in  1947  [7].  For  the  case  where  only  noisy  bearings  are  known 
of  a  single  stationary  emitter,  the  setup  as  described  by  Stansfield  is  loosely  pictured  in 
Figure  SI 

Typically,  the  problem  to  be  solved  is  that  while  the  bearing  lines  due  to  two  receivers 
intersect  with  no  ambiguity,  in  general  the  lines  of  three  receivers  do  not  intersect.  Instead, 
they  form  a  so-called  “cocked  hat” ,  a  triangle  within  which  the  emitter  might  be  thought 
probably  to  lie.  In  fact,  in  such  a  situation  where  the  bearing  line  drawn  is  the  centre 
of  some  probability  distribution,  there  is  only  a  1  in  4  chance  that  the  emitter  will  lie 
within  the  cocked  hat  triangle,  despite  the  fact  that  somewhere  within  the  triangle  seems 
to  be  the  best  place  to  estimate  its  position  in  the  absence  of  any  other  information. 
The  best  way  out  of  this  difficulty  is  to  put  any  method  of  locating  the  emitter  on  a 
good  statistical  footing.  We  assume  the  bearings  are  spread  according  to  some  known 
probability  distribution  in  bearing  angle,  peaked  around  the  true  bearing;  and  we  require 
the  most  probable  emitter  position  given  this  set  of  measured  bearings.  This  is  known  as 
the  maximum  a  posteriori  estimate  of  the  emitter’s  position.  So  we  need  to  calculate  the 
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Figure  3:  Stansfield  scenario 


probability  that  the  emitter  is  at  some  position  (x,y)  given  that  bearing  set,  and  vary  x 
and  y  over  the  whole  region  where  the  emitter  could  possibly  be,  until  we  find  a  maximum. 

It  proves  easier,  but  just  as  useful,  to  follow  a  slightly  different  approach  known  as 
the  method  of  maximum  likelihood.  In  this  we  take  the  true  position  of  the  emitter  to 
be  the  one  most  likely  to  have  given  rise  to  the  measured  bearings.  That  is,  we  calculate 
the  probability  that  the  measured  bearings  will  be  what  they  are,  for  all  possible  emitter 
positions,  and  choose  the  position  maximising  this  probability. 

That  is,  the  maximum  a  posteriori  approach  maximises 

proh{emitter  at  (x,y)  \  hearing  data) ,  (3-7) 

while  the  maximum  likelihood  approach  maximises 

proh{bearing  data  \  emitter  at  {x,y)) .  (3-8) 

These  two  probabilities  are  related  via  Bayes’  rule. 


Stansheld’s  original  calculation  is  based  around  the  following  maximum  likelihood  anal¬ 
ysis.  Begin  by  referring  to  Figure  Q]  on  page  [51  The  measured  bearing  is  just  the  exact 
bearing  plus  some  noise: 

Pk  =  Ok{x,y)  +  Vk  ■  (3.9) 

Provided  the  bearing  errors  are  normally  distributed,  the  maximum  likelihood  probability 
will  be  a  product  of  gaussian  functions: 

p(/3i . .  ./3n|6>i . .  .6>n)  PC  expy~]  2^^  ■  (3.10) 

k 


Maximising  this  probability  is  equivalent  to  minimising  twice  the  sum  of  squares  in  the 
exponent,  this  being  known  as  the  “cost  function”  of  the  problem: 


ML  cost  function  = 


{Pk  -  Okf 


(3.11) 
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Stansfield  expressed  (3k  —  9k  as  a  function  of  the  emitter  position  {x,  y),  and  then  equated 
each  of  the  x,  y  partial  derivatives  with  zero.  This  can  be  done  in  different  ways,  depending 
on  what  level  of  accuracy  is  chosen  and  what  assumptions  are  made.  Different  assumptions 
will  lead  to  slightly  different  formulations  and  solutions  of  the  problem. 

While  Stansheld’s  original  paper  used  the  maximum  likelihood  approach,  it  was  writ¬ 
ten  very  obscurely.  It  was  rewritten  a  decade  later  by  Ancker  [8]  using  a  slightly  more 
transparent  approach.  The  calculation  is  complicated  by  the  fact  that  the  probability 
distributions  are  expressed  fundamentally  in  terms  of  angles,  so  that  when  converted  to 
distances  using  cartesian  coordinates,  the  distance  of  the  emitter  from  each  receiver  enters 
the  problem  in  a  mathematically  complicated  way.  This  is  not  a  problem  in  itself;  it  just 
makes  the  resulting  equations  harder  to  solve. 

Both  Stansheld  and  Ancker  tackled  this  difficulty  by  assuming  that  the  distance  from 
each  receiver  to  the  trial  position  of  the  emitter  does  not  change  as  this  trial  position  is 
varied.  This  approximation  is  viable,  because  the  required  probability  will  be  nowhere 
near  its  maximum  when  this  distance  has  begun  to  extend  to  regions  where  we  know  the 
emitter  is  certainly  not  present;  thus  we  are  really  only  considering  the  area  around  the 
peak  of  the  probability — and  this  is  precisely  the  only  region  in  which  we  were  interested 
anyway. 

The  so-called  Stansfield  solution  for  geolocation,  as  restated  by  Ancker,  writes  the 
emitter’s  best-estimate  position  in  an  expression  that  still  involves  these  unknown  dis¬ 
tances,  and  as  such  is  an  iterative  technique  using  a  batch  of  bearings.  Ancker  suggested 
that  initial  estimates  of  these  distances  be  taken  from  the  centre  of  the  largest  polygon 
formed  by  the  bearing  lines.  But  calculating  a  good  centre  for  the  largest  polygon  formed 
becomes  a  problem  in  itself  not  unlike  the  original.  Ancker  does  add  that  in  the  absence 
of  such  an  estimate,  we  might  resort  to  a  trial  and  error  approach. 

In  the  following  pages,  we  first  consider  two  alternative  formulations  of  Stansfield’s 
problem,  each  having  its  own  method  of  solution:  the  Cartesian  Pseudo-Linear  Estima¬ 
tor  (CPLE),  and  the  Gauss-Newton  technique.  These  are  “batch  routines”,  in  that  they 
process  a  set  of  points  to  arrive  at  some  estimate  of  the  emitter  position.  This  is  followed 
by  discussion  and  examples  of  a  more  modern,  recursive,  approach  to  geolocation,  where 
only  the  latest  bearing  measurement  is  used  to  update  the  estimated  emitter  position. 
Sorenson  [9]  gives  further  background  to  the  use  of  least  squares  in  a  DE  context. 


3.6  A  Least  Squares  Primer 


The  discussions  on  the  CPLE  and  Gauss-Newton  routines  that  follow  all  use  the  formalism 
of  the  theory  of  least  squares,  which  we  review  here.  Suppose  we  have  taken  n  measure¬ 
ments  Zi,  which  in  a  noiseless  world  would  fit  linearly  to  a  set  of  m  parameters  xi, ,  Xm, 
using  well-defined  coefficients  that  comprise  an  n  x  m  matrix  H: 


Xl 

=  H 

Xm 

or  just 

z  =  Hx . 


(3.12) 


(3.13) 
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In  general,  we  have  an  abundance  of  measurements:  n  >  m  so  that  H  is  taller  than  square. 
In  a  noisy  world  the  elements  Xi  are  required  to  be  found  through  the  inexact  relation 

z  =  Hx  + noise,  (3-14) 

where  the  noise  cannot  be  known  exactly.  In  general  this  set  of  equations  is  overdeter¬ 
mined,  but  we  can  find  a  best  fit  for  a;  in  a  least-squares  sense,  by  choosing  x  to  minimise 
the  distance  between  the  two  points  in  n-dimensional  space,  Hx  and  z.  This  is  equivalent 
to  choosing  the  x  that  minimises  \Hx  —  zp  (called  a  “cost  function”,  which  equals  the 
squared  length  of  the  n-dimensional  noise  vector).  This  is  done  by  setting  the  gradient  of 
the  cost  function  to  zeroU 

V{|ifa;-z|2}  =  V[{Hx- z)\Hx- z)] 

=  V  [x^H^Hx  -  2z^Hx  +  z*z] 

=  2x^H^H  -2z^H  .  (3.16) 

Equating  the  last  expression  with  zero  defines  the  least-squares  solution  x: 

X  =  {H^H)-^H^z  =  H*z  ,  (3.17) 

where  is  the  pseudo  inverse  of  H.  Numerical  packages  such  as  Matlab  do  not  just 
calculate  by  inverting  and  transposing  as  in  (3.17).  Numerical  matrix  inversion  can 
be  hazardous  for  near-singular  matrices,  so  more  sophisticated  methods  are  used,  such 
as  singular  value  decomposition.  Matlab’s  pinv  function  makes  use  of  singular  value 
decomposition  (pinv(H)  =  77^). 


Error  Analysis 


In  practice,  each  measurement  might  be  subject  to  some  error  cj^.  Assuming  gaussian 
statistics,  it  is  then  more  meaningful  to  minimise  not  the  length  of  77a;  —  z,  but  rather  to 
weight  each  of  its  components  by  the  corresponding  inverse  variance.  That  is,  we  should 
minimise 

“Ik 


E 


(Hx-z)l 


k  k 

This  is  identical  to  minimising  a  new  cost  function; 

[Hx  -  zfp-^{Hx  -  z) 

with  respect  to  x,  where  P  is  a  diagonal  matrix  of  variances: 

O' 


p  = 


0 


CfZ, 


Use  the  following  easily-verified  results: 

Va;  {ad x)  =  ad  and  Va,  (a:*Tx)  =  x^{A  -\-  A) . 


(3.18) 


(3.19) 


(3.20) 


(3.15) 
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The  calculation  is  almost  identical  to  that  of  (j3.16|).  with  the  following  result: 

x  =  (3.21) 

If  all  the  errors  are  the  same,  i.e.  cJi  =  •  •  •  =  then  the  variance  matrix  P  cancels  from 
this  expression,  leaving  us  with  the  usual  pseudo  inverse  expression  (|3.17|). 

3.6.1  Total  Least  Squares:  When  the  Errors  are  More  Complicated 

On  a  more  advanced  note,  the  theory  of  how  best  to  manipulate  matrices,  when  real-world 
noisy  numbers  are  involved,  can  be  used  to  improve  upon  the  least-squares  analysis  of  the 
preceding  pages.  The  relevant  theory  is  known  as  singular  value  decomposition,  and  deals 
with  efficient  ways  to  express  matrices  in  terms  of  factors  that  are  robust  in  numerical 
work.  Robustness  can  be  necessary  when  small  numbers  are  involved  in  the  calculations, 
that  lie  at  the  edge  of  the  computing  machine’s  inbuilt  precision.  Such  numbers  might 
be  created  in  a  numerical  algorithm  being  used;  matrix  inversion  is  particularly  prone  to 
this.  Geolocation  using  singular  value  decomposition  is  discussed  in  Appendix  [Bl 


3.7  Cartesian  Pseudo-Linear  Estimator  Approach  (CPLE) 

Our  first  variant  of  the  maximum  likelihood  approach  (as  is  shown  later)  relies  on  the 
use  of  the  orthogonality  of  vectors  specifying  bearings,  as  shown  in  Figure  @1  In  essence 
what  the  method  does  is  locate  a  best  point  of  intersection  of  the  bearing  lines  with  a 
least-squares  fit.  The  calculation  is  simple  enough  to  be  done  nonrecursively,  and  does 
not  require  an  initial  estimate  of  the  emitter  position.  CPLE  is  in  essence  equivalent  to 
the  Stansfield  algorithm,  but  the  following  description  of  it  is  conceptually  simpler,  and 
unlike  the  usual  exposition  of  Stansfield,  is  easily  extended  to  moving  emitters. 

The  great  utility  of  the  CPLE  is  that  although  a  batch  process — a  once-only  processing 
of  a  given  set  of  data — it  is  still  accurate  enough  to  compete  with  the  other  routines 
discussed  ahead.  It  can  be  run  on  an  initial  data  set  and  then  used  again  on  sets  of 
subsequent  data,  but  it  can  also  be  used  to  provide  an  initial  estimate  of  where  the  emitter 
is,  in  order  to  supply  one  of  the  following  algorithms  with  their  required  initial  estimates 
of  the  emitter’s  position  and  state  of  motion. 

In  this  section  the  CPLE  method  is  outlined  and  the  moving-emitter  case  is  covered,  in 
which  the  emitter  is  modelled  as  having  at  most  a  constant  acceleration  in  two  dimensions. 
Extension  to  three  dimensions  and  higher  derivatives  of  position  is  straightforward. 

Suppose  that  the  emitter  position  is  a  vector  s{t).  There  are  n  receiver  positions,  where 
the  position  is  a  vector  Write  all  vectors  as  columns  and  work  in  cartesian  coor¬ 

dinates.  The  emitter  has  initial  position  Sg,  initial  velocity  Vq  and  constant  acceleration  a 
(all  vectors),  so  that 

s{t)  =  Sq  +  v^t  +  ^at‘^  .  (3.22) 

Group  the  three  constant  vectors  that  we  wish  to  find  into  one  column  vector  x.  This  has 
six  elements  if  we  are  confined  to  the  plane,  since  each  of  the  vectors  SQ,VQ,a  then  has 
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North 


Figure  4-  Cartesian  Pseudo-Linear  Estimator  approach 


two  elements: 


X  = 


a 


(3.23) 


In  tandem  with  x,  introduce  a  matrix  A{t)  that  serves  to  absorb  the  emitter  dynamics: 


A{t)^ 


1  0  t  0  %  0 
0  1  0  t  0  § 


so  that  s{t)  =  A{t)  X  . 


(3.24) 


This  leaves  the  state  x  of  initial  conditions  remaining  to  be  found. 

The  specifying  of  x  has  recast  the  problem  into  a  stationary  viewpoint.  At  some 
time  t,  receiver  k  makes  a  bearing  measurement,  and  notes  that  the  direction  from  it  to 
the  emitter  is  given  by  some  not-quite- known  noisy  vector  bi^{t),  plus  some  noise  i^i^{t).  (In 
practice  the  lengths  of  bj^{t),  i^i^{t)  are  completely  unknown  or  not  even  very  well-defined, 
but  that  is  of  no  consequence  to  the  calculation.)  We  must  estimate  x  given  the  set  of  all 
the  and  6^.  For  each  k  we  can  write 


s{t)  =  rk{t)  +  bf^{t)  +  i^f^{t) . 


(3.25) 


Now,  every  6^  has  two  unit  vectors  that  are  orthogonal  to  it.  These  turn  out  to  be 
useful,  although  since  we  are  working  in  two  dimensions  for  simplicity,  only  one  of  these 
will  be  needed:  call  it  5^,  so  that  b^  ■  bf,  =  0.  The  utility  of  b^  is  that  although  we 
have  no  knowledge  of  the  length  of  6^.,  we  do  know  its  direction,  which  is  enough  to 
specify  bj:  exactly.  That  being  the  case,  projecting  all  vectors  along  the  bj:  direction  will 
then  eliminate  the  unknown  .  This  is  effected  by  forming  the  dot  product  of  both  sides 
of  (3.25)  with  b^,  giving 


bi:  ■  s  =  bu  ■  Vi,  -\-  bu  ■  Vu  . 

fv  rv  i\j  t\>  i\j 


(3.26) 
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If  the  receiver  takes  its  bearing  measurement  at  time  (which  times  need  not  all  be 
different),  then  the  last  equation  becomes 


hi^A{t^)x  =  ht  -Tk  +  ht  -  Vk- 


(3.27) 


Writing  this  for  each  k  gives 

b^^A{tn)_ 

This  will  match  the  basic  equation  (|3.14|).  provided  we  identify 


bi  ri 


•  V 

.  n  '  n. 


X  +  residual  noise  term. 


bj;^A{ti) 

'bi  ■  rP 

H  = 

,  Z  = 

b~^  •  T 

L  n  '  nj 

in  which  case  the  least-squares  solution  can  be  written  down  immediately: 


X 


=  H*z  = 


'bi*A(tiy 

# 

bi  ■  ri 

bi^A(tn)_ 

pi  ■  rn_ 

(3.28) 


(3.29) 


(3.30) 


Because  the  matrix  H  is  built  from  noisy  bearings,  the  method  of  total  least  squares  might  be 
useful  here;  so  rather  than  apply  (1B10|  or  (13.301,  we  might  wish  to  try  (BIT).  In  fact,  it  turns  out 
that  total  least  squares  gives  10-20%  smaller  CEPs  than  the  ordinary  least-squares  approach. 


The  CPLE  is  a  simple  locating  method,  yet  gives  very  good  results,  even  when  using  just 
the  minimum  number  of  data  points  required  (for  low  noise).  This  number  can  be  small. 
Consider,  for  example,  the  stationary  emitter  case  in  two  dimensions,  where  just  two 
bearings  are  needed  (since  these  specify  two  lines  with  a  well-defined  intersection  where 
the  emitter  can  lie).  Note  that  two  bearing  lines  are  also  sufficient  in  three  dimensions, 
because  again  we  are  only  finding  their  intersection;  since  each  line  needs  two  numbers 
to  specify  it  (altitude  and  azimuth),  we  have  four  numbers  in  total — more  than  enough 
information  to  solve  for  the  three  emitter  coordinates.  Alternatively,  in  three  dimensions 
we  are  using  two  unit  normals  for  each  bearing  line,  so  that  with  just  two  lines  (13.26)  is 
really  four  equations,  and  hence  s  is  already  overspecified. 

Being  a  batch  process,  the  CPLE  routine  can  use  all  of  the  data  to  produce  a  one-off 
estimate  of  the  emitter  position,  and  incoming  new  bearings  cannot  be  processed  individu¬ 
ally  to  refine  this  estimate.  But  in  practice  the  CPLE  gives  a  good  estimate  with  very  few 
data  points,  and  so  can  always  be  run  on  a  subset  of  data  points  that  includes  the  latest 
one.  Also,  its  simplicity  means  it  can  be  used  to  generate  an  estimate  quickly — which  can 
then  be  used  to  seed  one  of  the  methods  described  in  the  next  sections.  Unlike  the  CPLE, 
these  methods  all  require  an  initial  estimate  of  the  emitter  state. 


16 


DSTO-RR-0319 


Error  Analysis 

Because  the  CPLE  is  a  straightforward  application  of  the  least-squares  technique,  the  way 
to  calculate  the  error  in  the  estimated  emitter  position  (xe,ye)  is  known  from  standard 
statistical  theory.  If  the  bearing  error  is  a  constant  cr,  then  the  relevant  covariance  matrix 
is  [3110]: 

cov(xe,?/e)  =  (3.31) 

As  discussed  in  Appendix  0,  if  we  calculate  this  matrix  and  use  it  to  draw  an  error 
ellipse  centred  on  the  estimate  point,  the  result  is  usually  correlated  well  with  the  actual 
performance  of  the  algorithm  in  terms  of  where  it  estimates  the  emitter  to  be  on  repeated 
runs.  This  error  ellipse  also  tends  to  correlate  well  with  the  Cramer-Rao  lower  bound 
ellipse — but  not  always.  It  might  be  much  larger  than  the  Cramer-Rao  ellipse  but  can 
also  be  somewhat  smaller.  It  does  tend  to  have  about  the  same  size  and  orientation. 
Presumably  the  departure  from  Cramer-Rao  is  caused  by  a  bias  in  the  CPLE. 


3.7.1  Relating  the  CPLE  to  the  Method  of  Maximum  Likelihood 


The  least-squares  solution  to  the  CPLE  turns  out  to  be  a  simplified  form  of  the  method 
of  Maximum  Likelihood.  We  can  show  this  by  comparing  the  CPLE  cost  function  to  the 
maximum  likelihood  cost  function.  Note  from  (3.271 13.281 13.29)  that  \Hx  —  z\'^  is  the  sum 
of  squares  of  the  noise  terms  Writing  as  the  true  distance  from  receiver  k  to 

the  possibly- moving  emitter  at  time  this  sum  of  squares  becomes  a  sum  of  squares  of 
perpendicular  distances: 

CPLE  cost  function  =  \Hx  —  z|^  =  ^ 

k 

k 

^  (3.32) 

k 


In  contrast,  the  maximum  likelihood  cost  function  was  calculated  on  page  TT  to  be: 


ML  cost  function  = 


(/3fc  -  (^k 


crt 


(3.33) 


This  compares  well  with  the  CPLE  cost  function,  provided  that  (a)  our  lack  of  knowledge 
of  the  distances  means  we  are  prepared  to  set  all  of  these  to  be  constant,  so  that  they 
factor  out  of  (3.32):  and  (b)  the  bearing  errors  are  all  the  same,  so  can  be  factored  out 
of  (3.33).  The  CPLE  does  indeed  stand  on  a  very  solid  statistical  footing. 


3.8  Gauss— Newton/MLE  Algorithm 

Given  a  set  of  bearings,  and  wishing  to  use  a  least-squares  ht  to  the  set  of  parameters 
required,  a  very  straightforward  and  widely-known  approach  again  applies  the  maximum 
likelihood  idea,  but  this  time  to  the  angular  difference  between  true  and  observed  emit¬ 
ter  positions.  It  has  attracted  the  name  maximum  likelihood  estimation,  and  the  usual 
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solution  method  used  is  the  Gauss-Newton  method  of  following  a  gradient  to  locate  a 
minimum  point.  These  are  generic  names,  but  both  have  been  attached  to  this  particu¬ 
lar  method,  which  is  termed  the  Gauss-Newton  algorithm  in  this  report  and  elsewhere. 
Further  discussion  can  be  found  in  the  reference  by  Foy  }llj.  There,  Foy  compares  the 
Gauss-Newton  and  Kalman  routines  (described  later),  noting  that  Gauss-Newton  is  more 
desirable  due  to  its  simplicity.  However,  like  the  GPLE,  the  Gauss-Newton  algorithm  is  a 
batch  technique,  and  so  hts  into  a  different  class  of  routines  than  does  the  Kalman  hlter. 

Begin  by  relating  measured  to  exact  bearings: 


h  —  ^ki^)  +  ■ 


(3.34) 


We  wish  to  recast  this  equation  in  the  matrix  language  of  Section  13.61  so  as  to  solve  it 
in  a  least-squares  sense.  One  way  to  do  this  first  estimates  the  emitter  state  (3.23)  to  be 
some  Xq,  and  then  Taylor  expands  (3.34)  around  this  estimate  to  hrst  order: 


where  the  gradient  is 


Pk 

-  ^fc(®o) 

+  V0,( 

Xq){x  - 

®o)  +  ^k  ! 

(3.35) 

dh 

dOk 

dh 

dh 

dOk  1 

(3.36) 

^(«oJ 

d(soy) 

^Ky).  ■ 

With  n  measurements  taken,  write  (3.35)  as 


- 1 

1 

-  ^ 

O 

_ 1 

'v9,{x,y 

1 

o 

e 

•  • 

1 

_ 1 

y9nixQ)_ 

{x-Xq)  +  V>, 


(3.37) 


to  be  solved  for  x.  Writing 


'v9yxy 

'Pi 

O 

1 

H  = 

,  z  = 

y9n{xQ)_ 

Pn 

o 

1 

(3.38) 


we  see  that  this  equation  now  matches  (3.14).  so  that  we  can  immediately  write  down  an 
updated  state  estimate  xi  with  a  least-squares  approach: 


Xl  =  Xn  + 


Wi(*o)- 

# 

1 

1 

-  ^ 

O 

_ 1 

_V6»„(®o)_ 

- 1 

o 

e 

•  • 

1 

_ 1 

(3.39) 


This,  the  Gauss-Newton  algorithm,  processes  a  batch  of  measurements  and  iterates  to 
converge  (hopefully)  to  the  actual  emitter  state.  Simulations  show  that  the  method  of 
total  least  squares  does  not  work  anywhere  near  as  well  for  the  Gauss-Newton  algorithm 
as  does  normal  least  squares. 

Use  the  notation  from  the  previous  pages,  and  set  to  be  the  emitter  position  at  the 
time  of  measurement  k.  (That  is,  we  are  free  to  model  the  emitter  as  having  e.g.  constant 
velocity,  constant  acceleration,  and  so  on.)  If  we  set 


Axi^ 

Aw 


k. 


=  Su  —  r 


k  5 


(3.40) 
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then  the  true  bearings  are  given  by 


(3.41) 


In  that  case,  if  the  emitter  motion  is  modelled  by  e.g.  constant  acceleration,  then  the 
gradient  with  respect  to  Sg^, . . .  ,ay  can  be  written  down: 


(3.42) 


where  12x2  = 


The  new  estimate  xi  will  lead  to  new  values  of  z  and  H,  so  that  we  must  iterate  until 
hopefully  the  estimates  converge  to  within  some  tolerance — although  there  is  no  guarantee 
that  this  limit  will  be  the  actual  emitter  position.  (Again,  a  bias  is  present  in  this  tech¬ 
nique,  just  as  occurs  in  the  CPLE.)  In  practice,  we  step  between  iterations  by  an  amount 
larger  than  that  given  in  (3.39):  that  is,  we  introduce  a  gain,  being  some  number  whose 
size  is  found  by  experiment  (but  is  roughly  1),  and  multiply  this  gain  by  the  step  in  (3.39). 
In  general,  increasing  the  gain  speeds  up  convergence — but  if  the  steps  are  too  large  then 
the  algorithm  will  become  oscillatory  or  unstable. 

In  writing  code  for  these  routines,  we  need  to  be  aware  that  the  lone  inverse  tangent  function 
needs  to  be  supplemented  with  a  constant  in  order  to  return  angles  within  the  full  360°  range. 

Of  course,  the  constant  differentiates  to  zero  so  need  not  be  explicitly  written  here — although  it 
certainly  must  be  inserted  into  code  that  implements  the  routine.  But  even  this  is  not  enough. 

In  the  actual  processing,  we  are  ultimately  comparing  measured  bearings  with  theoretical  values. 
Depending  on  the  convention  chosen,  there  must  always  be  a  critical  emitter  position  where  the 
theoretical  value  jumps  by  360°  as  the  emitter  is  moved  slightly.  This  lack  of  continuity  will  lead 
to  instability  in  the  routine,  unless  the  code  is  written  carefully  to  take  account  of  this  fact. 

Error  Analysis 

Although  some  sort  of  error  analysis  can  be  introduced  into  the  one-run  CPLE,  a  similar 
approach  to  the  iterative  Gauss-Newton  algorithm  is  not  clear-cut.  Each  iteration  is 
a  simple  least-squares  process;  however,  what  results  is  not  an  estimate  of  the  emitter 
position,  but  rather  an  estimate  of  the  step  size  needed  to  move  away  from  the  initial 
estimate,  together  with  the  error  in  this  step  size,  should  we  wish  to  calculate  this  error 
by  the  same  approach  as  was  used  for  the  CPLE.  But  we  cannot  just  start  out  with  some 
initial  estimate  with  its  own  error,  and  then  simply  add  more  errors  as  several  more  steps 
are  taken,  since  the  errors  are  presumably  not  independent:  this  approach  would  yield 
a  final  error  much  greater  than  what  is  actually  observed  through  running  simulations. 
Gauss-Newton  can  give  quite  accurate  results  even  with  a  very  inaccurate  initial  estimate, 
but  by  no  means  always — it  can  also  diverge  wildly  when  the  other  methods  described  here 
are  well-behaved. 

Spingarn  [3]  writes  the  error  for  Gauss-Newton  as  our  (3.31).  For  scenarios  in  which 
Gauss-Newton  does  converge,  the  hnal  error  is  described  fairly  well  by  the  ellipse  plotted 
from  the  covariance  matrix  of  (3.31) .  This  error  ellipse  also  usually  approximately  matches 
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the  Cramer-Rao  bounding  ellipse.  An  alternative  representative  error  for  Gauss-Newton 
is  often  found  by  computing  the  CRLB  at  the  final  estimate  of  the  state  x. 

Note  that  like  the  CPLE,  Gauss-Newton  is  also  a  batch  routine,  in  that  it  requires 
a  new  observation  to  be  added  to  at  least  some  of  the  previous  ones  before  running  the 
algorithm.  In  this  sense  it  might  be  considered  wasteful,  since  at  least  some  old  measure¬ 
ments  need  to  be  reprocessed.  The  need  for  a  streamlined  approach  led  historically  to  the 
Recursive  Least  Squares  class  of  algorithms. 


3.9  Recursive  Least  Squares  (RLS) 

Recursive  Least  Squares  algorithms  form  a  class  of  update  algorithms  that  process  the 
latest  data  point  only.  Such  techniques  have  been  known  for  the  last  half  century,  and  two 
of  the  more  widely  known  ones  are  considered  here.  The  Hrst  is  a  procedure  that  dates 
back  to  1950  if  not  earlier  [12],  and  this,  at  least  in  the  form  used  by  Godard  in  1974,  has 
come  to  be  called  simply  the  Recursive  Least  Squares  algorithm,  and  is  very  popular  [13] . 
The  second  algorithm  is  the  famous  Kalman  filter,  due  to  Kalman  in  1960  [14] . 

The  principle  behind  these  routines  is  that  since  we  have  already  made  use  of  all  except 
the  latest  data  point,  any  step  in  the  algorithm  that  uses  all  of  the  data  should  be  re- 
expressible  in  terms  of  past  values  and  the  latest  data.  This  results  in  an  update,  using 
only  the  latest  data,  of  the  last  state  estimate.  Because  of  this,  the  amount  of  data  needing 
to  be  manipulated  is  kept  at  a  constant  small  size:  e.g.  the  matrix  H  of  the  previous  two 
routines  is  now  always  composed  of  just  one  row,  without  growing  a  new  row  with  each 
incoming  data  point,  as  it  would  do  in  the  previous  routines. 

The  usual  RLS  routine  as  presented  here  is  derived  in  [15] ,  although  it  is  here  extended 
to  the  constant  acceleration  case  for  the  emitter.  As  usual,  the  extension  to  higher  order 
motions  is  straightforward.  The  matrices  H  and  z  of  the  previous  sections  are  now  replaced 
by  the  latest  measurements  only,  giving  them  sizes  1x6  and  1x1  respectively.  In  analogy 
to  (]3.38[)  we  write 

H  =  Ve{x),  z  =  (5-e{x).  (3.43) 

The  procedure  followed  using  RLS  is,  however,  more  complicated  than  that  of  Gauss- 
Newton.  It’s  written  in  the  following  way  so  as  to  be  suggestive  of  the  Kalman  filter, 
discussed  next: 

1.  Start  with  some  estimate  of  the  matrix  of  variances  P  [as  in  (j3.20j)].  equal  to  some 
large  number  h  times  the  6x6  identity.  A  good  choice  of  h  is  crucial  and  needs  to 
be  determined  empirically. 

2.  Ghoose  a  number  r  between  0  and  1,  again  to  be  determined  empirically.  This  dehnes 
a  sort  of  fading  memory,  in  the  sense  that  it  allows  more  weight  to  be  given  to  more 
recent  measurements. 

3.  Galculate  H,z  from  (3.43). 

4.  Galculate  a  gain  K  =  PH^/{r  +  HPH^). 
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5.  Now  update  the  emitter  position: 

X — >x  +  Kz.  (3.44) 

6.  Finally  update  P,  and  then  return  to  step  3: 

P^{hx6-KH)P/r.  (3.45) 

This  algorithm  enables  us  to  keep  updating  the  latest  estimate  of  the  emitter  position,  by 
only  ever  just  incorporating  the  latest  measurement  as  it  arrives. 


3.10  Kalman  Filter 

The  Recursive  Least  Squares  algorithm  is  not  always  stable,  and  several  attempts  have 
been  made  to  improve  on  it.  The  best  known  is  due  to  Kalman  (14],  and  for  linear  gaussian 
models  the  Kalman  filter  is  the  only  one  in  wide  use  [16].  Its  basic  principles  are  explained 
well  by  Sorenson  [9],  whose  notation  is  used  here.  The  algorithm  has  been  rederived 
recently  in  an  accessible  way  by  Challa  and  Koks  [17],  using  the  ideas  of  Bayes  theory. 
The  algorithm  is  conventionally  called  the  “Kalman  filter” ,  because,  like  Recursive  Least 
Squares,  it  filters  the  data  to  give  a  predicted  estimate  of  the  state  at  each  timestep.  This 
differs  from  methods  such  as  CPLE  and  Gauss-Newton,  which  are  smoothers,  in  that  they 
operate  on  a  set  of  data  to  produce  just  one  state,  which  is  then  used  to  predict  the  future 
motion  of  the  emitter. 

A  basic  premise  underlying  Kalman  theory  is  that  the  system  being  observed  evolves 
in  some  way  able  to  be  modelled.  In  the  case  of  tracking,  the  emitter  might  be  modelled 
with  e.g.  constant  velocity,  so  that  its  state  comprises  its  position  and  velocity.  In  general 
its  state  at  the  time  of  the  measurement  can  be  written  as  some  vector  x^  whose 
evolution  we  model  as  follows: 

Xk+i  =  Fk^k  +  Gk^k  +  U). .  (3.46) 

Here,  Ff,  is  some  evolution  matrix  allowing  us  to  specify  a  model  of  how  the  emitter’s 
state  changes,  is  the  “process”  noise,  a  term  that  models  how  the  emitter’s  motion 
might  depart  from  what  our  model  assumes  through  the  term.  For  the  purpose  of 
investigation  we  have  built  the  filter  by  assuming  the  emitter  moves  in  some  complicated 
way,  by  putting  the  effects  of  acceleration  (and  its  derivatives,  if  need  be)  into  the  evolution 
matrix  F^,  as  is  done  in  [I8l  |T9].  We  can  also  model  any  jitter  in  the  target  by  another 
kind  of  acceleration  term,  but  rather  than  explicitly  putting  it  into  we  include  it  in  a 
stochastic  way  through  the  matrix  that  multiplies  the  noise.  Finally,  is  any  extra 
non-noise  term  required,  such  as  a  shift  in  coordinates. 

The  state  Xj^  is  not  observed  directly;  rather,  we  make  a  measurement  z^,: 

^k  =  ^k^k  F'^kPUk-  (3.47) 

Here,  describes  the  measurement  process  in  terms  of  the  system  state,  w^,  is  any  noise 
introduced  by  the  measurement,  and  is  any  extra  term  required  as  part  of  the  model. 
The  process  and  measurement  noises  are  assumed  white  and  gaussian: 

~  Af(0,  Qfc) ,  ~  iV(0,  Rfc) .  (3.48) 
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The  standard  terminology  also  defines: 

=  predicted  value  (“predictor”)  of  determined  after  measurement  k  —  1 
=  best  estimate  (“filtered  estimate”)  of  determined  after  measurement  k. 

(3.49) 

The  Kalman  filter  also  allows  us  to  predict  the  error  in  £Cfc|fc_i  and  via  the  covariance 
matrices: 

Pk\k-i  =  ~  ®fc|A:-i)(same)*}  =  covariance  of  error  in  predicted  estimate  of  Xf, 

=  E{(a;^  —  $fc|fc)(same)*}  =  covariance  of  error  in  filtered  estimate  of  Xj^. 

(3.50) 

The  Kalman  filter  is  implemented  through  the  following  procedure: 


1.  Estimate  £o[o>  -folo^ 

For  k  =  1  onwards,  calculate  the  following: 

2-  Pk\k-i  =  Pk-i\k-i  Fk-i  +  G'fc-i  Qk-i  G'fc-i 

3-  +  '^k-i 

4.  K,  =  Hi  [H,  Hi  +  R,)-^ 

5-  ®A:|A:  =  ^k\k-l  +  ^k\k-l  ~  l/fc) 

6-  Pk\k  =  {I  —  ^k^k)  Pk\k-1  (3.51) 

V _ ^ _ J 


Note  that  and  can  be  calculated  offline  if  the  process  and  measurement 

matrices  are  known  for  all  time,  as  they  often  are;  and  this  makes  the  linear  Kalman  filter 
very  fast  to  implement. 

The  Kalman  filter  does  not  have  the  same  sort  of  gain  as  found  in  Gauss-Newton,  nor 
a  fading  memory  parameter  as  found  in  RLS.  It  simply  combines  the  state  inferred  from 
the  latest  measurement,  with  the  current  state  of  the  system  as  predicted  from  the  other 
measurements  to  date,  using  appropriate  weightings;  and  it  does  this  in  a  statistically 
robust  and  optimal  way. 

Despite  its  many  components,  there  is  no  common  notation  in  general  use  for  the 
Kalman  filter.  For  example,  the  noise  vector  and  scalar  v  and  w  are  commonly  called 
w  and  V. 

Extended  Kalman  Filter  (EKF) 
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The  linear  Kalman  filter  gives  the  best  possible  fit  for  linear  problems  with  gaussian  noise, 
i.e.  where  the  state  parameters  and  the  measurements  taken  at  time  k  +  1  are  linear 
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combinations  of  the  state  parameters  at  time  k.  But  no  such  claim  can  be  made  for  the 
nonlinear  case,  such  as  is  encountered  in  geolocation  through  (j3.41j).  The  modification 
made  to  the  linear  filter  to  incorporate  nonlinearity  gives  rise  to  the  Extended  Kalman 
filter  (EKF).  We  consider  here  two  types:  the  first  uses  cartesian  coordinates,  while  the 
second,  both  more  effective  and  more  complex,  uses  a  modified  set  of  polar  coordinates. 

To  extend  the  linear  filter’s  application  to  these  nonlinear  cases,  start  with  analogues 
of  the  linear  equations  (13.461 13.471) : 

=  fki^k)  +  9ki'^k) 

^k  =  KiXk)+'Wk,  (3.52) 

and  linearise  these  assuming  Xf,  —  to  be  small,  using  expressions  like 

fki^k)  -  fki^kik)  +  Fk  ■  (Xk  -  Xk\k) 

=  -ffc  •  (3.53) 

Ff^  is  a  matrix  whose  row  is  the  gradient  of  the  component  of  /^,  with  respect  to  the 
set  of  state  parameters.  Similar  equations  hold  for  and  and  define  (It  seems 

to  be  normal  to  expand  both  and  around  x^^  but  around  this  is  just 

an  arbitrary  convention,  but  we  have  followed  suit.)  The  filter  is  applied  in  the  following 
way: 


1.  Estimate  SqIO)  -fo|0)  Qki  before. 

Then  for  k  =  1  onwards: 

2-  Pk\k-1  =  Pk-l\k-l  Pk-l  +  F^k-l  Qk-l  F^k-l 

3-  Xk\k-1  =  /fe-l(®fc-l|fc-l) 

5-  ®A:|A:  =  +  ^k  [^k  ~  ^A:(®fc|A:-l)] 

6-  Pk\k  —  {I  ~  P-k^k)  Pk\k-l  (3.54) 


Two  problems  arise  in  a  nonlinear  scenario.  The  first  is  that  since  we  are  always 
expanding  around  the  latest  estimate  of  the  state  vector,  the  matrices  P^j^_i,  Kj^,  and 
Pfcjfc  can  no  longer  be  calculated  offline,  and  this  adds  greatly  to  the  time  taken  to  run 
the  filter.  The  second  problem  is  that  since  the  covariance  matrix  P^j^  depends  on  77^ — 
which  is  now  being  calculated  from  the  latest  estimates  of  the  state  vector  in  a  similar 
way  to  (3.53) — the  covariances  inherit  the  inaccuracy  of  this  state  vector,  and  so  can  no 
longer  be  relied  upon  to  provide  the  true  errors  in  the  results.  Hence  there  appears  to  be 
no  good  method  of  estimating  the  accuracy  of  the  result,  apart  from  measuring  the  spread 
of  estimated  emitter  positions  or  tracks  based  on  a  Monte  Carlo  approach.  Worse,  the 
inaccuracies  in  i-^j^  now  feed  back  into  the  Kalman  equations  and  can  produce  instabilities. 
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Cartesian  EKF 


The  tracking  problem  can  be  formulated  in  the  following  way  for  the  case  of  an  emitter  with 
constant  acceleration  in  cartesian  coordinates.  Remember  that  in  the  RLS  and  Kalman 
approaches,  we  don’t  characterise  the  state  of  the  emitter  by  the  constants  parametrising 
its  motion  in  the  way  that  we  did  for  the  CPLE  and  Gauss-Newton.  Rather,  the  emitter’s 
state  vector  Xf,  reflects  its  state  at  the  instant  k.  If  we  write  its  position  at  k  as 
(hoping  that  the  use  of  both  and  does  not  cause  confusion!),  then  the  state  vector 
is  (with  dots  denoting  time  rates  of  change,  and  a  transpose  used  to  save  vertical  space): 


Xi 


T  =  [xk  Vk  Xk  Vk  Xk  Vk] 
Next  we  write  the  state  equation 

'1  0  At  0 

0  1  0  At 

0  0  1  0 

0  0  0  1 

0  0  0  0 

0  0  0  0 

Writing 


^k+l  — 


AtV2 

0 

At 

0 

1 

0 


0 

AtV2 

0 

At 

0 

1 


=  Fk^k 


'^k 

.^Vk. 

Vk. 

—  r 


k  > 


then  the  measurement  model  is  just  the  bearing  seen  by  the  receiver: 

Zu  =  tan“^  - — -  +  noise, 

^  Axu 


=  Kalman  h^{x^) 


in  which  case  we  have 
hki^k)  -  K{xk\k-i)  + 


dx. 


dhk 

dVk 


dK 

dx. 


dh,  dh,  dh. 


dy.  dx,  dy_ 


k. 


(3.55) 


(3.56) 


(3.57) 


(3.58) 


{Xk-Xk\k-i),  (3.59) 


=  Kalman  H, 


where 

0  0  0  0]-  (3-60) 

The  algorithm  can  now  be  implemented,  and  this  has  been  done  using  Matlab  as  described 
ahead  in  Appendix  [Cl 


Challa  and  Faruqi  [20]  have  made  a  comparison  of  two  approaches  to  using  this  filter 
for  a  stationary  emitter.  They  use  two  versions  of  the  state  vector  The  first  is  the 
one  used  here,  i.e.  the  emitter  position,  a  2  x  1  vector.  The  other  version  they  use  is  a 
4x1  vector  composed  of  the  emitter’s  position  and  its  velocity  relative  to  the  receiver, 
with  all  components  expressed  as  fractions  of  the  range.  They  then  show,  by  running 
several  simulations  of  Spingarn’s  example  [3],  that  the  second  choice  of  state  vector  leads 
to  better  results:  faster  convergence,  and  also  more  accurate  convergence  in  the  presence 
of  noise.  This  begs  the  question  of  what  actually  is  the  best  choice  of  state  vector,  and 
this  is  a  large  area  of  ongoing  research  in  the  geolocation  community,  that  has  not  been 
investigated  in  this  report. 
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Error  Analysis 

Like  Gauss-Newton,  errors  in  the  Recursive  Least  Squares  and  Kalman  routines  are  also 
difficult,  if  not  impossible,  to  evaluate.  When  the  Kalman  filter  is  applied  to  linear  tasks, 
the  errors  in  its  state  estimates  are  given  by  the  entries  of  the  covariance  matrix  Pk\k]  but 
for  a  linear  approximation  to  the  nonlinear  problem  of  geolocation  as  used  here,  the  Pk\k 
matrix  loses  some  of  its  covariance  meaning  even  though  its  role  in  the  Kalman  algorithm 
is  unchanged.  See  the  discussion  on  page  [231 

Since  both  RLS  and  the  Kalman  filter  are  iterated  as  each  new  data  point  comes  in, 
the  covariance  evolves  over  time.  Thus,  while  e.g.  the  error  in  the  x  component  of  position 
as  seen  in  P^^k  might  get  smaller  to  an  acceptable  final  value,  the  y  position  error  may  well 
start  out  by  getting  smaller,  but  might  soon  begin  to  grow;  it  might  even  become  much 
larger  than  the  final  error  in  the  estimated  placement  of  the  emitter.  Or,  it  might  shrink 
to  zero.  As  in  the  Gauss-Newton  case  and  to  some  extent  the  CPLE  also,  it  could  well  be 
that  the  best  way  to  calculate  a  state  estimate  error  is  through  a  Monte  Carlo  approach, 
by  processing  repeated  runs  of  the  receiver  flight  path  and  looking  at  the  scattering  of 
estimates  to  the  emitter  position  that  result. 


Modified  Polar  Coordinate  EKF 


The  Cartesian  EKF  is  not  always  stable  and  can  give  biased  estimates,  as  shown  explicitly 
in  [20] .  Variations  have  been  proposed  that  are  based  instead  on  polar  coordinates  [201  [21] . 
We  consider  one  such  well  known  way  here:  the  Modified  Polar  coordinate  EKF  algorithm 
(MPEKF)  [I9l.  \22[  \21\.  The  increased  stability  provided  by  this  algorithm  is  paid  for 
initially  by  a  more  complicated  choice  of  state  vector.  In  the  notation  of  Arulampalam  [19], 
if  f3j^  is  the  bearing  angle  of  the  emitter  as  seen  by  the  receiver  measured  clockwise  from 
north  (i.e.  Arulampalam’s  =  7r/2  —  our  (5^  in  Figure  Ql),  and  if  Arulampalam’s  is 
the  range  to  the  emitter  measured  from  the  receiver,  then  the  state  vector  corresponding 
to  (I3.55D  is  set  to  be 


Xk 


(3.61) 


Unlike  the  previous  state  vectors  described,  this  is  independent  of  whether  or  not  the 
emitter  accelerates,  since  we  are  now  modelling  the  effects  of  acceleration  only  by  including 
a  nonzero  matrix  the  covariance  of  the  process  noise.  So  this  form  for  the  state 
vector  suffices  no  matter  how  the  emitter’s  motion  is  modelled.  The  more  complex 
entries  of  the  state  vector  mean  that  the  implementation  of  the  filter  is  more  involved: 
refer  to  [19]  for  the  details  of  calculating  the  initial  matrices  that  model  covariances  and 
how  the  system  evolves. 

As  shown  in  Appendix  [C]  the  performance  of  the  MPEKF  is  dependent  on  a  good 
initialisation,  which  of  course  is  something  that  cannot  necessarily  be  achieved  in  practice. 
The  filter’s  performance  has  also  been  analysed  in  [2l],  and  is  shown  there  to  be  stable 
and  asymptotically  unbiased. 
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4  Review  of  Selected  Papers  in  the  Field 

Stansfield’s  original  geolocation  paper  became  much  cited  in  the  literature,  although 
the  main  effect  it  had  was  to  spark  interest  in  finding  efficient  ways  to  solve  the  geolocation 
problem.  At  the  time,  nomograms  and  specially  graded  least-squares  rulers  were  used  for 
working  with  charts  overdrawn  with  bearing  lines,  but  the  advent  of  more  computing 
power  also  saw  the  rise  of  more  complex  methods,  such  as  the  Kalman  filter,  that  depend 
for  their  use  on  computers. 

The  Kalman  filter  is  often  seen  as  fundamental  to  tracking  (and  for  example  is  used  in 
Spingarn’s  widely-quoted  paper  [3]),  but  the  simpler  Gauss-Newton  method  is  still  quite 
comparable  in  its  performance.  This  is  especially  true  with  more  powerful  computers  that 
allow  the  fast  reprocessing  of  data  as  per  Gauss-Newton’s  needs.  For  example,  Poirot 
and  M'^Williams  }23]  apply  Gauss-Newton  successfully  to  bearings  as  measured  on  Earth 
(i.e.  using  latitude  and  longitude).  Rao  and  Reddy  [24]  also  deal  with  Earth  using  a  least- 
squares  approach  that  makes  use  of  new  combinations  of  the  various  bearing  parameters, 
and  by  combining  this  new  routine  with  a  Kalman  filter,  have  produced  a  more  efficient 
algorithm.  This  highlights  an  important  point  in  the  search  for  more  efficient  algorithms: 
by  solving  the  problem  using  new  and  unusual  combinations  of  the  conventionally  used 
variables,  an  increase  in  efficiency  can  sometimes  be  obtained  (as  is  done  in  [20]);  but  an 
approach  to  producing  these  unusual  combinations  is  not  clear. 

Elliott  and  others  [16]  have  recently  produced  an  algorithm  that  generalises  the  Kalman 
filter  and  appears  to  perform  much  better,  but  their  analysis  is  somewhat  obscure,  and 
the  algorithm  is  not  considered  in  this  report. 

Gavish  and  Weiss  [25]  compare  the  performance  of  a  least-squares  algorithm  with 
the  more  involved  Stansfield  algorithm.  They  state  a  well-known  theorem  that  says  the 
Maximum  Likelihood  Estimator  is  unbiased,  and  can  achieve  the  Gramer-Rao  lower  bound 
in  its  accuracy  provided  the  number  of  measurements  is  large  enough.  By  running  scenarios 
using  both  a  least-squares  approach  and  the  Stansfield  algorithm,  they  find  that  the  least- 
squares  bias  vanishes  as  the  number  of  measurements  increases,  while  the  Stansfield  results 
are  biased  and  this  bias  does  not  vanish,  and  can  even  grow,  as  the  number  of  measurements 
increases. 

Mahapatra  takes  a  different  approach  to  geolocation  in  general  [26] .  He  suggests  flying 
a  receiver  flight  path  such  that  the  bearing  of  the  stationary  emitter  is  held  constant;  the 
parameters  of  the  flight  path  then  give  the  emitter  location  without  actually  using  the  DF 
bearing.  In  this  way,  the  moving  receiver  is  free  to  concentrate  on  keeping  an  accurate  hx 
on  the  emitter  without  needing  to  record  exactly  what  this  fix  is.  The  down  side  to  this  is 
that  a  very  special  path  geometry  needs  to  be  flown,  which  of  course  does  not  necessarily 
suit  the  goals  of  a  given  flight  mission. 

The  effect  of  bearing  bias  in  passive  geolocation  has  also  been  studied  by  Gavish  and 
Fogel  [27].  who  consider  a  constant  angle  added  to  the  emitter  bearing.  As  we  have  done 
in  Section  13.41  they  define  the  notion  of  observability  to  mean  that  in  the  absence  of  noise 
but  perhaps  with  this  bias,  the  emitter’s  position  is  still  able  to  be  uniquely  determined. 
They  show  that  if  there  is  any  possibility  of  an  additive  bearing  bias  being  present,  then 
the  receiver  should  not  fly  a  circle  that  has  any  chance  of  passing  though  the  emitter’s 
position. 
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Their  argument  runs  as  follows.  Suppose  that  to  the  right  hand  side  of  (j3.41j)  we 
allow  a  constant  bearing  bias  to  be  present,  and  disregard  the  noise  so  that  we  are  only 
considering  the  systematic  error  of  the  added  bias.  Then  the  notion  of  observability  means 
there  is  only  one  emitter  position  allowed  as  a  solution  to  this  equation.  But  suppose  that 
on  the  contrary  there  are  two  distinct  solutions,  perhaps  with  different  biases.  Then  it’s 
straightforward  to  show  that  the  track  that  gives  these  solutions  is  a  circle  passing  through 
the  emitter’s  actual  position.  So  the  moral  is  that  if  we  want  to  be  assured  of  not  having 
any  ambiguity,  then  we  should  not  invite  difficulties  by  flying  any  such  circle  in  the  first 
place. 

Gavish  and  Fogel  also  consider  the  Cramer-Rao  lower  bound  and  how  it  is  affected 
by  a  bias.  In  particular,  they  consider  the  case  of  a  constant  velocity  receiver  flying 
approximately  broadside  to  an  emitter.  The  bias  added  to  the  bearings  is  no  longer  a 
constant,  but  rather  gaussian-distributed  around  zero  mean  with  some  specified  standard 
deviation.  They  calculate  a  Cramer-Rao  ellipse  for  a  number  of  different  track  lengths, 
assuming  a  2°  standard  deviation  noise  and  both  with  and  without  a  3°  s.d.  bias.  What 
they  hnd  is  that,  as  might  be  expected,  the  error  ellipse  is  larger  for  the  biased  case,  and 
as  a  simple  rule  of  thumb  it  could  be  said  that  typically  the  addition  of  the  bias  means 
that  the  major  axis  of  the  no-bias  ellipse  becomes  the  minor  axis  of  the  biased  ellipse,  with 
the  shape  of  the  ellipse  being  roughly  unchanged. 

Reports  of  Healy  and  of  Beasley  &  Miles.  Here  we  describe  two  longer  reports 
that  give  an  overview  of  stationary-emitter  geolocation  in  general,  as  well  as  comparing 
different  methods. 

(1)  Healy’s  thesis.  The  first  is  a  thesis  by  Healy  [28].  In  this  he  compares  the 
Gauss-Newton  and  RLS  methods  for  the  case  of  both  stationary  and  moving  emitters,  with 
his  main  effort  going  toward  quantifying  the  algorithms’  effectiveness.  After  discussing  the 
algorithms  and  the  general  problem  to  be  solved,  he  runs  each  for  a  variety  of  scenarios 
and  calculates  several  things:  the  speed  of  convergence,  misadjustment  accuracy  (how 
well  the  algorithms  pinpoint  the  emitter),  computational  complexity  (number  of  floating 
point  operations  required  for  the  calculations),  and  robustness  to  initial  data  and  different 
conditions  of  data  generation. 

A  selection  of  Healy’s  results  have  been  verihed  or  compared  using  the  Matlab  pro¬ 
gramme  described  in  our  Appendix  O  One  of  Healy’s  initial  results  is  that  the  speeds  of 
convergence  of  G-N  (Gauss-Newton)  and  RLS  are  similar  once  the  RLS  has  overcome  its 
slow  start,  during  which  time  it  needs  to  stabilise  the  matrix  P  in  (|3.45).  In  fact  this  is, 
however,  quite  dependent  on  how  we  implement  the  G-N  algorithm.  Healy  divides  the 
data  points  into  batches  of  ten,  overlapping  by  five,  and  runs  the  algorithm  separately  on 
each  batch.  In  the  Matlab  programme  of  Appendix  (Cl  this  has  not  been  done;  instead  the 
algorithm  is  always  run  on  the  entire  set  of  points  since  we  have  enough  computational 
speed  to  justify  this.  What  this  means  is  that  just  one  calculation  is  done  at  the  end  of  the 
flight  path.  Doing  this  and  recording  the  number  of  floating  point  operations  required  by 
Matlab  to  implement  the  various  algorithms,  it  turns  out  that  for  fairly  straightforward 
scenarios,  Gauss-Newton  takes  about  three  times  as  many  floating  point  operations  as  do 
both  RLS  and  Kalman.  Healy’s  results  also  agree  with  this. 
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On  his  page  77,  Healy  states  that  G-N’s  misadjustment  accuracy  is  fairly  well  inde¬ 
pendent  of  its  gain.  By  running  our  Matlab  programme  we  find  that  although  G-N  does 
generally  converge  for  each  of  the  two  gains  that  Healy  uses  (0.3  and  1),  it  might  oscil¬ 
late  wildly  and  iterate  a  very  large  number  of  times  before  approaching  a  limit  value.  In 
agreement  with  Healy,  we  find  that  RLS  is  very  sensitive  to  gain  changes,  as  discussed  in 
Appendix  \0^ 

Healy  also  finds  that  RLS  and  G-N  are  about  equally  robust  when  applied  to  various 
situations  of  moving  emitters.  On  his  page  87,  Healy  writes  that  the  G-N  algorithm’s 
stability  is  not  affected  by  the  choice  of  gain.  Although  it  is  certainly  true  that  changing 
the  gain  will  not  make  the  algorithm  diverge,  a  change  may  well  cause  a  convergence  to  a 
very  different  final  value  if  the  gain  is  too  small,  or  cause  it  to  oscillate  heavily  about  the 
true  value  if  the  gain  is  too  large — and  also  the  numbers  of  iterations  required  for  even 
these  scenarios  are  heavily  dependent  on  the  gain  choice.  So  it  would  be  naive  to  infer 
that  just  because  an  algorithm  is  stable,  its  results  can  be  trusted. 

Many  of  Healy’s  trials  use  a  moving  emitter,  for  which  he  has  success  with  the  RLS 
algorithm.  Understandably  his  Gauss-Newton  routine  fails  here,  being  only  designed  for 
the  stationary  case.  In  Section  13.8  we  have  tailored  the  G-N  routine  to  cope  with  moving 
emitters. 


(2)  Report  by  Beasley  and  Miles.  The  second  report  available  is  by  Beasley  and 
Miles  [2] .  This  provides  a  general  overview  of  both  direction  finding  and  the  associated 
statistical  processing  for  geolocation,  such  as  maximum  likelihood.  It  then  focuses  on 
TDOA,  using  least-squares  analysis  with  gaussian  errors. 

The  report’s  modelling  section  describes  the  calculation  of  error  ellipses,  that  are  asso¬ 
ciated  with  geolocation  from  a  constant  velocity  aircraft  employing  a  small  TDOA  receiver 
at  each  wingtip.  As  expected,  these  ellipses  show  that  errors  are  minimised  for  broadside 
emitters  (the  closer  the  better),  with  the  down-range  error  increasing  rapidly  as  the  emitter 
angle  moves  away  from  broadside.  This  is  a  reasonable  result,  since  the  time  differences 
of  signal  reception  are  then  becoming  less  sensitive  to  distance.  In  that  case  too,  the 
corresponding  cross  range  errors  also  increase,  but  are  much  smaller  in  size.  Perhaps  the 
main  results  from  the  modelling  section  in  [2]  are  that  the  down-  and  cross-range  errors 
can  be  approximated  as  follows: 


Down  range  oc  ■  ; 

LdVNB^TP 


Cross  range  oc  —  , 

dVNB^TP 


(4.1) 


where 


R  =  range  to  emitter 
d  =  baseline  of  sensor  pair 
B  =  emitter  bandwidth 
P  =  emitter’s  effective  radiated  power. 


L  =  receiver  flight  path  length 
N  =  number  of  bearings  taken 
T  =  length  of  signal  wave  train 


(Although  these  expressions  are  specifically  for  a  TDOA  scenario  with  two  receivers,  they 
are  referred  to  on  our  page  06]  in  connection  with  a  more  general  proposed  rule  of  thumb.) 
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Also  included  in  Beasley  and  Miles’  report  are  studies  of  ways  to  beat  systematic  errors. 
One  approach  for  estimating  any  bias  uses  reference  emitters  at  known  locations.  Beasley 
and  Miles  do  this,  but  for  their  example  the  method  offers  only  a  small  improvement  in 
accuracy.  Another  suggested  way  of  cancelling  bias  is  by  the  use  of  race  track  paths  for 
the  receiver,  in  which  it  flies  along  a  straight  line  and  then  back  again,  with  the  aim  of 
cancelling  any  errors  introduced  by  the  emitter’s  being  always  on  the  same  side  of  the 
plane. 

The  report  also  discusses  the  physical  characteristics  required  of  the  receiver  antenna 
for  good  reception  of  the  signal,  and  briefly  discusses  the  relevant  signal  processing  re¬ 
quirements.  It  finishes  by  describing  DF  work  carried  out  by  the  United  Kingdom  Defence 
Evaluation  and  Research  Agency  (DERA)  in  Malvern,  Defford,  Portsdown,  and  Reading. 


One  area  not  well  described  in  the  geolocation  literature  is  the  shape  of  tracks  flown, 
not  just  by  emitters  but  also  by  receivers.  Spingarn  [3]  and  Challa  &  Faruqi  [201  129] 
use  a  constant  velocity  receiver  with  stationary  emitter;  Arulampalam  [19]  and  Lindgren 
(fe  Gong  [30]  use  a  zig-zagging  receiver  flying  at  constant  velocity  on  each  leg,  following 
a  constant  velocity  emitter.  Farina  [22]  and  Galkowski  &  Islam  [31]  do  much  the  same 
but  with  a  single  change  in  receiver  direction.  This  simplicity,  especially  in  older  litera¬ 
ture,  is  presumably  the  result  of  the  authors’  not  having  fast  enough  computers  at  their 
disposal  to  be  programmed  easily  for  more  complex  motion.  It  is  addressed  by  the  Mat- 
lab  programme  described  in  Appendix  El  which  incorporates  receivers  and  emitters  flying 
arbitrary  tracks — at  least  insofar  as  these  are  describable  using  a  series  of  waypoints. 


Moving  Emitters.  Much  less  space  in  the  literature  has  been  devoted  to  the  more  in¬ 
teresting  subject  of  moving  emitters,  although  several  authors  have  discussed  geolocating 
a  constant  velocity  emitter.  Lindgren  and  Gong  [30]  describe  a  Kalman  filter  with  a  simple 
receiver  motion  consisting  of  just  two  constant  velocity  legs,  and  find  good  convergence; 
although  they  concede  that  precisely  how  the  receiver  should  move  to  maximise  the  al¬ 
gorithm’s  performance  is  a  matter  for  further  study.  Levine  and  Marino  [32]  discuss  a 
similar  problem,  proving  theorems  about  observability  and  the  advantages  of  knowing  the 
emitter’s  distance,  but  they  carry  out  no  simulations  using  actual  data. 


Other  Kalman  Filters.  Arulampalam  has  published  a  report  [19]  in  which  he  compares 
the  performance  of  various  Kalman  filters.  Besides  the  Cartesian  and  Modified  Polar 
EKFs,  he  also  considers  others:  the  Pseudo-Linear  Estimator  EKF,  Modified  Gain  EKF, 
and  range-parametrised  versions  of  the  Cartesian  and  Modified  Polar  EKFs. 

The  first  two  of  these  other  filters  are  built  on  modifications  to  or  in  Section  |3.10i 
The  range-parametrised  versions  of  the  CEKF /MPEKF  are  built  essentially  by  running 
the  usual  CEKF /MPEKF  Liters  in  parallel  for  different  range  initialisations,  and  then 
combining  the  results  in  a  way  that  uses  bayesian  analysis  to  provide  weightings.  This  can 
all  be  computationally  intensive,  but  in  practice  some  of  the  parallel  calculations  have  a 
negligible  enough  weighting  as  to  exclude  them  from  needing  to  be  performed  at  all. 


Arulampalam’s  work  shows  the  following.  The  Pseudo-Linear  Estimator  is  generally 
only  as  good  as,  or  worse  than,  the  others;  although  along  with  the  MPEKE  it  outperforms 
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the  others  in  giving  a  heading  error  that  quickly  drops  to  zero  with  time.  On  the  other 
hand  it  can  give  a  poor  range  error.  It  is  not  considered  further  in  this  report.  The 
performance  of  the  Modified  Gain  EKF  is  generally  similar  to  that  of  the  CEKF  and  also 
is  not  considered  further  here.  The  two  range-parametrised  trackers  are  found  to  be  of 
most  efficacy  only  when  there  is  little  or  no  knowledge  of  the  initial  range.  However,  in 
Arulampalam’s  scenarios  they  are  perhaps  50%  better  at  estimating  emitter  speed  than 
their  CEKF /MPEKF  counterparts,  while  their  range,  azimuth,  and  heading  errors  differ 
little  from  those  of  the  usual  CEKF /MPEKF.  These  range-parametrised  algorithms  have 
also  been  excluded  from  further  study  in  this  report. 


Receiver  Position  Errors.  Little  work  has  been  published  on  the  case  where  the  re¬ 
ceiver  positions  are  themselves  subject  to  error.  Ancker  [8],  in  his  rewrite  of  the  Stansfield 
algorithm,  considered  the  case  where  the  linear  errors  in  the  receiver  positions  are  gaus- 
sian  and  small  compared  to  the  emitter  distances,  but  large  enough  to  be  comparable  with 
linear  errors  corresponding  to  the  bearing  error.  He  calculated  an  effective  bearing  error 
by  simply  adding  the  two  linear  errors  from  the  receiver  position  and  the  emitter  location. 
This  requires  knowledge  of  the  emitter  range  (as  discussed  previously),  which  in  practice 
can  be  estimated  initially  and  then  again  at  each  step  of  an  iterative  solution.  Wax  [33] 
considers  the  same  problem  and  generalises  it,  although  the  analysis  is  somewhat  obscure. 

It  is  more  realistic  for  the  receiver  positions  not  to  have  normally  distributed  errors,  but 
rather  to  be  subject  to  some  offset.  Receiver  positions  that  are  calculated  using  GPS  have 
a  very  small  error.  More  important  is  the  gyro  error  of  the  inertial  navigation  system,  since 
this  cannot  be  calibrated  using  GPS.  It  can  be  modelled  in  a  first  analysis  as  increasing 
linearly,  although  it  is  also  subject  to  an  oscillation.  But  this  error  is  also  very  small:  a 
typical  value  for  our  two-dimensional  scenario  would  be  a  0.001°  hr“^  increase  in  yaw,  so 
that  this  would  manifest  as  a  bias  of  the  same  amount  to  each  bearing  measured.  We  have 
not  considered  such  small  values  any  further. 


5  Geolocation  For  a  Triangle  of  Receivers 

In  this  section  we  consider  a  more  restricted  problem  than  previously.  Given  a  geometry 
consisting  of  three  stationary  receivers  placed  in  an  equilateral  triangle,  using  bearings- 
only  geolocation  for  a  stationary  emitter,  how  good  a  geolocation  can  be  done? 

Our  start  point  is  a  calculation  of  the  Cramer-Rao  bound,  being  the  theoretical  best- 
obtainable  accuracy.  The  calculation  is  done  here  in  detail  for  its  pedagogical  value.  Next 
we  run  several  simulations  of  various  geolocation  algorithms,  and  measure  parameters  such 
as  the  CEP.  These  parameters  are  then  plotted  against  the  emitter’s  distance  from  the 
triangle  along  the  direction  in  which  the  Cramer-Rao  results  indicate  that  we  can  expect 
the  worst  estimates.  Figure^ shows  the  basic  geometry  always  used. 


5.1  Calculating  a  Cramer-Rao  Set  of  Bounds 
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A  set  of  error  ellipses  that  describe  the  Cramer-Rao  lower  bound  for  a  multiple  receiver 
setup  is  produced  in  this  section.  As  discussed  in  Section  13.21  the  underlying  principle  is 
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Figure  5:  Calculation  of  the  Cramir-Rao  bound  for  three  receivers  observing 
one  emitter 


that  the  inverse  of  the  Fisher  Information  matrix  is  the  best  covariance  matrix  obtainable 
for  an  unbiased  estimator  of  the  emitter  location. 

The  Fisher  Information  matrix  is  produced  in  the  following  way.  Suppose  we  have  a 
set  of  n  stationary  receivers,  each  observing  a  stationary  emitter  (as  in  Figure  where 
n  =  3).  The  exact  bearing  at  receiver  i  is  denoted  6*: 


bi  =  tan 


ye  -  yi 


+  quadrant-dependent  constant. 


(5.1) 


We  define  the  actual  (i.e.  noisy)  bearing  observed  by  receiver  i  as  hi  with  error  e*  =  bi  —  bi. 
These  errors  are  gaussian-distributed  around  zero  with  standard  deviation  dj: 


P{^i)  = 


1 


cr, 


exp 


Z£! 

2a} 


(5.2) 


The  procedure  for  calculating  the  Fisher  Information  matrix  begins  with  the  probability 
density  for  all  of  the  receivers.  They  are  independent,  so  we  can  write 


p(ei,...,en)  =  p{£i) . .  .p{en)  ■ 

As  in  (3.3).  calculate  the  negative  log-likelihood  function  defined  by 


(5.3) 


L(a;e,  Ve)  =  -  lnp(ei, . . . ,  en)  =  ^  In  (aiV^j  + 


2ai 


2  • 


(5.4) 


The  Fisher  Information  matrix  itself  is  then  the  expected  value  of  a  matrix  of  partial 
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derivatives  of  this  negative  log-likelihood: 


■  d^L 
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(5.5) 


The  expectations  are  taken  with  respect  to  the  gaussian  in  (|5.2[) — which  is  of  course  an 
even  function,  in  which  case  the  following  holds: 


j 

dbj  dbj 
dxe  dye 


dbi  dbi 
dxe  dye 

[dye 


(5.6) 


The  partial  derivatives  are  found  from  fl5.1j).  producing 

{Ve-Vif  -{Xe  -  Xi){ye  -  yi) 

J  -  ~  -  Vi)  {Xe  -  Xif 


crf[{xe  -  Xi)^  +  {ye  -  yi)^]‘ 


(5.7) 


The  Cramer~Rao  theory  states  that  is  the  best  covariance  matrix  obtainable  from 
an  unbiased  estimator  of  {xe,ye)-  The  c-sigma  uncertainty  ellipse  around  each  emitter 
position  is  then  given  by  the  following  plot: 


[x  -  Xe 


y-ye]i 


X  —  Xe 

y- ye 


C  . 


(5.8) 


If  we  take  the  scenario  in  Figure  ^  with  1°  bearings  errors  (i.e.  all  CTj  =  1°),  then  a  plot 
of  3-sigma  ellipses  is  shown  in  Figure  [61  We  see  the  expected  threefold  symmetry,  as 
well  as  the  fact  that  the  best  results  (smallest  ellipses)  occur  when  the  overall  angular 
spread  of  the  receivers  is  largest,  when  seen  from  the  emitter  (which  lies  at  the  centre 
of  any  given  ellipse).  Thus  any  emitter  at  the  centre  of  the  northernmost  ellipses  (and 
their  counterparts  at  bearings  of  120°  and  240°)  see  the  largest  receiver-subtended  angle, 
corresponding  to  these  ellipses  having  the  smallest  major  axes.  The  reverse  is  true  for  the 
southernmost  ellipses  and  their  counterparts.  Minor  axes  comparisons  are  not  so  easy  to 
make. 
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What  we  now  wish  to  do  is  plot  a  representative  dimension  of,  say,  each  la  Cramer- 
Rao  ellipse  for  the  worst  case  scenario  (where  two  receivers  are  equally  further  distant  than 
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Figure  6:  Cramer-Rao  3a  ellipses  for  the  setup  of  Figure  [5i  The  receivers 
are  blue  crosses,  with  the  centre  of  each  ellipse  being  the  actual  position  of 
the  emitter  used  in  the  calculation  of  that  Cramer-Rao  ellipse.  We  have  used 
3a  instead  of  the  more  usual  la  purely  to  make  the  ellipses  large  enough  for 
easy  viewing 

the  closest  receiver  to  the  emitter),  as  a  function  of  the  emitter  distance  from  the  centre 
of  the  triangle.  This  worst  case  scenario  corresponds  to  an  emitter  being  at  a  bearing 
of  60°  in  Figure  [61  What  we  will  plot  is  the  circular  error  probable  (CEP),  being  the 
radius  of  the  circle  centred  on  the  true  emitter  position,  within  which  the  emitter  would 
be  estimated  to  lie  50%  of  the  time.  This  has  been  found  by  calculating  the  radius  of  the 
circle,  centred  at  the  emitter,  that  makes  a  domain  such  that  when  the  double  gaussian  of 
the  Cramer-Rao  ellipse  is  integrated  over  this  domain,  the  answer  is  0.5.  The  resulting 
plots  are  in  Figure  [7]  for  bearing  errors  of  1°  and  8°,  along  with  baselines  of  20  and  40  km. 


5.2  Long-Baseline  TDOA 

In  this  section,  we  wish  to  establish  an  accuracy  for  geolocation  using  the  long-baseline  time 
difference  of  arrival  technique,  again  for  the  above  equilateral  triangle  of  three  receivers 
drawn  in  Figure  El  Each  of  these  receivers  is  located  at  some  distance  R  from  the  origin, 
and  each  has  an  inherent  timing  error  of  r — meaning  that  the  times  they  register  are 
in  error  by  an  amount  drawn  from  a  normal  distribution  with  zero  mean  and  standard 
deviation  r. 

Suppose  that  the  signal  arrives  at  receiver  1  at  time  t,  at  receiver  2  at  time  t  -|-  Ati2, 
and  at  receiver  3  at  time  t  -|-  Ati^.  Then  the  radio  waves  have  travelled  at  speed  c 
for  an  extra  distance  cAti2  to  get  to  receiver  2  as  compared  with  receiver  1,  and  sim¬ 
ilarly  for  receiver  3.  So  if  the  unknown  emitter  location  is  {x,y),  with  the  receivers  at 
positions  {xi,yi)  — >  (x3,y3),  then  we  can  write,  for  di  — >  d^  being  the  distances  of  the 
receivers  from  the  emitter: 

di  =  i/ {xi  -  x)2  -h  {yi  -  yY  for  i  =  1  ^  3 
d2  =  di  -I-  cAti2 

d's  =  di  +  cAti3  .  (5.9) 
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Figure  7;  50%  CEP  for  Cramer-Rao  estimates,  applied  to  the  setup  of 
Figure  [5^  with  the  emitter  at  a  60°  bearing  as  seen  from  the  centre  of  the 
receiver  triangle.  Top  row:  hearing  error  1° ;  baselines  20km  (left)  and  fOkm 
(right).  Bottom  row:  bearing  error  8°;  baselines  20km  (left)  and  fOkm 
(right) 


Given  time  differences  Ati2,Ati3,  equation  (15. 9f)  can  be  solved  numerically  to  find  the 
emitter  position.  See  the  box  on  the  facing  page  for  how  to  do  this. 

Conversely,  if  we  start  with  knowledge  of  the  emitter  position,  then  this  TDOA  can 
be  used  to  establish  a  geolocation  accuracy  for  this  technique  when  the  timing  error  r  is 
present:  by  simulating  the  geolocation  many  times  and  building  up  a  spread  of  estimated 
positions  for  the  emitter.  First  we  place  the  emitter  at  some  well  defined  position,  and 
calculate  values  for  Ati2,Ati3.  Thus  each  geolocation  simulation  can  be  made  by  first 
perturbing  these  values,  by  adding  to  each  of  them  a  random  number  chosen  from  the 
normal  distribution  A^(0,r^);  and  then  treating  the  new  “noisy”  time  differences  as  simu¬ 
lated  data  from  which  we  calculate  an  estimated  position  for  the  emitter  by  solving  fl5.9j). 
Doing  this  many  times  gives  a  spread  of  estimated  emitter  positions,  and  we  can  calculate 
a  CEP  for  this  spread,  relative  to  the  true  position  of  the  emitter. 
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Solving  Nonlinear  Simultaneous  Equations 

How  do  we  solve  the  set  of  equations  (15. 9p?  Write  it  as 


d2{x,y)  -  di{x,y)  -  cAti2  =  0 

d3{x,y)  -  di{x,y)  -  cAtis  =  0  (5.10) 


There  are  two  equations  in  x  and  y  here.  With  x  =  [x  y]*,  write  the  first 
equation  as  fi{x)  =  0  and  the  second  as  f2(x)  =  0.  For  some  Xq  not  too 
far  from  x,  we  can  Taylor-expand  them  both  to  first  order,  writing  (and 
remembering  that  V/i  and  V/2  are  row  vectors): 


hix) 

'h{xo) 

1 

V/i(xo)  (x  -  Xo) 

f2ix) 

M^o) 

"T 

_V/2(xo)  (x  -  Xo) 

= 

fiixo) 

-t 

■V/i(xo)l  .  _  X 
[v/2(xo)J 

=  J{xo) 


If  fi{x)  =  f2{x)  =  0  by  definition,  then 


fiixo) 

f2ixo)^ 


+  J{xo)  {x 


Xq)  ~  0, 


(5.12) 


or 


X~Xo-J  ^(xo) 


fl{xo) 

Mxo)^ 


(5.13) 


This  equation  can  be  used  iteratively:  the  initial  estimate  of  the  emit¬ 
ter  position  is  Xq,  and  the  new  estimate  is  given  by  the  right  hand  side 
of  (M). 


When  this  is  done  using  (j5.13j),  we  can  get  a  feel  for  how  the  geolocation  is  a  function 
of  the  baseline  and  timing  accuracies.  The  timing  accuracy  is  more  crucial,  by  which  is 
meant  that  if  a  given  long  baseline  TDOA  system  has  a  certain  geolocation  accuracy,  and 
we  decide  to  halve  the  baseline,  then  in  order  to  keep  the  stated  geolocation  accuracy, 
we  must  do  more  than  halve  the  number  specifying  the  timing  accuracy.  This  is  seen  in 
the  setup  of  Figure  ^  (page  I3l]).  with  a  triangle  side  of  10km  and  an  emitter  100  km  to 
the  east.  The  receivers  have  a  Ins  timing  accuracy.  Following  the  procedure  outlined  in 
the  previous  paragraph,  the  50%  CEP  is  calculated  to  be  2.3  km.  Now,  if  we  reduce  the 
triangle  side  length  by  a  factor  of  ten  to  1  km,  does  the  accuracy  required  to  keep  the 
same  CEP  now  also  reduce  by  a  factor  of  10  to  become  100  ps?  No:  the  required  timing 
accuracy  turns  out  to  be  about  1  ps.  So  in  order  to  reduce  the  size  of  the  TDOA  array, 
we  need  a  drastically  improved  timing  capability.  This  is  shown  in  Figure  [81 

Note  that  for  TDOA,  the  geolocation  accuracy  is  quite  geometry-dependent.  For  the 
above  case  of  a  10km  triangle  side  with  Ins  timing,  if  the  emitter  is  100km  north  of  the 
triangle  (see  Figure  [9),  then  the  2.3  km  CEP  reduces  to  just  0.1km.  In  that  case,  scaling 
the  triangle  down  by  a  factor  of  10  places  less  of  a  demand  on  the  timing  accuracy  to  keep 
the  0.1km  CEP:  instead  of  Ips,  we  need  just  10  ps. 
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10  km  side 


1  ns  timing 


2.3  km  CEP 


1 00  km  away  . 


(i)  Original  scenario 


•  •  1  km  side 


•  1  ps  timing 


2.3  km  CEP 


100  km  away 


(ii)  Scaled-down  version 


Figure  8:  A  schematic  of  a  Monte  Carlo  set  of  estimated  emitter  positions. 
Beginning  with  a  10  km  triangle  and  1  ns  timing  (i),  the  50%  CEP  is  2.3  km. 

To  maintain  this  CEP  while  scaling  the  receiver  geometry  down  by  a  factor 
of  10  means  inereasing  the  timing  accuracy  by  a  factor  of  1000,  as  seen  in  (ii) 

Linear  Analysis 

Although  a  linear  analysis  of  TDOA  might  be  possible,  it  is  certainly  fraught  with  dif¬ 
ficulty.  The  reason  is  that  the  emitter  is  being  located  at  the  intersection  of  (at  least) 
two  hyperbolm,  and  any  linear  assumption  might  affect  the  position  of  the  distant  arms  of 
these  that  are  doing  the  intersecting.  Although  it  might  be  thought  that  as  the  emitter’s 
distance  increases,  a  linear  assumption  will  be  increasingly  accurate,  in  practice  a  compet¬ 
ing  effect  is  present:  the  distance  to  the  intersection  of  the  appropriate  hyperbolas  arms 
will  become  ever  more  sensitive  to  the  position  of  those  arms  as  they  become  more  and 
more  parallel.  So  even  a  slight  error  for  a  distant  emitter  can  produce  a  huge  error  in  the 
geolocation — even  in  the  absence  of  noise. 

A  good  example  of  this  occurs  in  the  following  way.  Suppose  we  first  assume  that  the 
incoming  rays  from  the  emitter  are  parallel,  and  then  use  TDOA  with  simple  trigonometry 
to  establish  a  direction  to  the  emitter.  This  method  gives  a  very  accurate  emitter  direction, 
but  it  cannot  be  used  in  practice  to  geolocate:  although  the  second  bearing  line  drawn 
with  the  aid  of  a  third  receiver  does  indeed  intersect  the  first,  the  placement  error  (for  no 
noise)  turns  out  to  be  prohibitively  large.  For  example,  in  the  scenario  of  Figure  [9)  the 
emitter  will  be  incorrectly  geolocated  at  about  twice  its  correct  distance  from  the  receiver 
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(i)  Original  scenario 


1  km  side 


•  10  ps  timing 


0.1  km  CEP 


100  km  away 


(ii)  Scaled-down  version 


Figure  9:  As  for  Figure  but  with  a  different  receiver  triangle  orienta¬ 
tion  (i).  The  50%  CEP  is  now  much  smaller:  0.1km.  To  maintain  it  now 
while  scaling  the  receiver  geometry  down  by  a  factor  of  10  means  we  need 
“only”  increase  the  timing  accuracy  by  a  factor  of  100  (ii) 


triangle.  In  the  final  analysis,  geolocation  with  TDOA  certainly  demands  good  baselines, 
timing,  and  numerical  accuracy. 


6  Results  of  Running  a  Matlab  Simulation 

This  section  describes  some  results  of  simulating  various  geolocation  scenarios.  These 
scenarios  were  used  as  a  testbed  for  comparing  the  various  geolocation  algorithms  described 
in  this  report.  The  algorithms  were  coded  into  Matlab  using  a  graphical  user  interface, 
or  “gui”,  designed  to  allow  reasonably  arbitrary  scenarios  to  be  set  up.  All  scenarios 
assumed  one  emitter  and  one  receiver,  with  either  moving  along  a  set  of  arbitrarily  laid- 
down  waypoints.  The  receiver  would  make  a  bearing  measurement  at  time  intervals  that 
could  also  be  specified  and  easily  changed.  These  measurements  were  subject  to  bearing 
noise  whose  error  could  be  input  and  changed  easily.  This  gui  is  described  in  Appendix  O 


The  primary  result  to  emerge  from  this  process  of  running  geolocation  simulations  is 
that  rules  of  thumb  for  efficient  geometries  appear  not  to  be  obtainable.  Simulations  can 
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always  give  a  feel  for  what  the  algorithms  are  capable  of,  and  what  their  foibles  are;  but 
simple  rules  of  thumb  have  not  arisen  from  any  of  the  work  described  in  this  section. 
Because  of  that,  the  software  is  mainly  useful  for  developing  a  feel  for  geolocation,  and 
is  available  from  the  author.  The  package  has  survived  a  Matlab  upgrade,  although  that 
upgrade  created  a  bug  that  was  only  fixed  with  difficulty.  Matlab  software  cannot  be 
assumed  to  be  backwards  compatible.  Also,  any  changes  caused  by  an  upgrade  will  not 
necessarily  cause  a  graceful  failure;  rather,  the  software  might  still  run,  but  give  slightly 
wrong  results — and  then  only  sometimes. 

The  output  from  the  following  command  is  a  simple  example  of  how  Matlab  can  give  an  empty- 
matrix  result  which  is  wrong  due  to  (necessary)  internal  rounding,  but  which  is  meant  to  be  “[3]”: 
find(  [0 . 5  :  0 . 1  :  0 . 8]  ==  0.7).  It  demonstrates  that  internal  rounding  errors  do  not  necessarily 
imply  a  result  that  is  incorrect  in  e.g.  the  tenth  decimal  place.  Rather,  they  can  be  completely 
wrong,  which  the  user  of  the  software  will  not  necessarily  be  aware  of.  Worse,  very  problematical 
to  the  code  used  in  the  gui  is  that  there  seems  to  be  a  new  Matlab  requirement  to  declare  double 
precision  in  places  not  listed  in  its  official  documentation,  and  which  was  not  the  case  in  previous 
versions;  and  the  explicit  use  of  double  precision  is  not  commonplace  in  Matlab  anyway.  The 
problem  is  that  not  declaring  double  precision  does  not  always  simply  lead  to  less  accurate  results; 
rather,  it  can  produce  meaningless  ASCII  characters  in  places  where  numbers  are  expected,  which 
leads  to  erratic  behaviour  in  the  gui.  (Submitting  modified  code  for  Matlab’s  find  command, 
along  with  a  request  asking  for  the  documentation  to  be  updated,  have  not  had  any  effect.) 

Some  care  needs  to  be  taken  with  coding  geolocation  routines  as  regards  the  trigonom¬ 
etry  used.  For  example,  Spingarn’s  well  known  work  [3]  uses  an  inverse  tangent  function; 
but  such  a  function  will  need  careful  coding  if  it  is  not  to  fail  in  some  situations.  An 
example  of  the  pitfall  can  be  set  up  with  two  receiver  trajectories  that  point  in  two  only 
slightly  differing  directions,  where  for  one  the  algorithm  succeeds,  while  for  the  other  it 
fails.  It  might  be  thought  that  since  (j3.41[)  uses  an  inverse  tangent,  then  all  we  need  do  is 
call  on  the  appropriate  one  of  Matlab’s  two  built-in  functions  for  this  purpose.  The  reason 
why  this  will  sometimes  fail  is  that  the  linearisation  of  the  routines  as  in  (j3.35p  expects 
the  z-vector  to  have  entries  that  become  smaller  with  each  iteration.  But  because  these 
entries  result  from  the  difference  of  two  angles,  we  might  have  a  situation  where  one  of 
the  angles  is  a  little  more  than  zero  (say  0.1  radians),  while  the  other,  given  via  an  inverse 
tangent  function,  is  a  little  under  27r:  say  6.2  radians.  While  these  two  angles  represent 
vectors  that  really  are  close  to  each  other,  their  numerical  difference  is  certainly  nowhere 
near  zero.  So  in  practice,  routines  for  dealing  with  angles  often  need  more  sophistication 
than  Matlab’s  inverse  tangent  functions  can  provide.  The  gui  software  has  been  written 
to  handle  this  situation. 

As  regards  computing  speed,  recent  increases  in  processor  speed  mean  that  the  3500 
bearing  computations  that  Healy  mentions  taking  15  minutes  for  RLS  in  1996-97  [28] 
now  take  just  a  few  seconds  on  a  standard  personal  computer.  This  is  a  good  indication 
that  geolocation  algorithms  can  now  be  more  realistically  used  than  they  could  even  a  few 
years  ago. 

In  an  effort  to  investigate  the  parameters  used  by  Healy  for  the  Gauss-Newton  and 
RLS  routines,  the  programme  used  in  this  report  was  first  run  using  a  fairly  straightfor¬ 
ward  initial  scenario:  a  50°  subtended  receiver  flight  path  as  seen  by  the  emitter,  emitter 
broadside  to  mid-path,  flight  path  length  100  units,  initial  estimate  of  emitter  position 
about  20  units  from  actual  position  and  somewhat  closer  to  the  flight  path  than  the  emit¬ 
ter.  A  very  convoluted  flight  path  with  multiple  loops  and  curves  was  also  used.  Bearing 
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error  was  set  at  1-2°,  which  is  a  very  typical  figure.  What  resulted  were  the  following 
parameter  choices: 

Gauss— Newton  gain  Healy  |28j  determines  the  best  value  by  asking  how  many  iter¬ 
ations  are  required  in  order  for  the  difference  in  successive  iterations  to  be  less 
than  0.001.  He  hnds  that  a  gain  of  1  gives  the  fastest  convergence,  although  appar¬ 
ently  has  not  considered  gains  higher  than  1.  When  running  the  Matlab  programme 
we  find  that  for  a  bearing  noise  of  1°,  the  number  of  iterations  required  to  get 
to  a  <  0.001  update  is  minimised  for  gain  values  in  the  range  0.95-1.05,  and  hence 
when  running  the  programme  the  gain  has  usually  been  set  equal  to  1. 

RLS  parameters  There  are  two  that  need  determining  (see  page  120]) : 

1.  The  constant  b.  The  value  of  this  is  not  well  discussed  in  the  literature,  but  is 
supposed  to  be  large  in  some  sense.  Healy  suggests  a  value  of  between  1  and 
10,  but  eventually  uses  both  150  and  0.1  with  a  moving  emitter.  The  value  of  b 
needs  to  be  set  in  conjunction  with  the  following  parameter  A. 

2.  The  value  of  A.  Healy  sets  this  equal  to  1.65k/{k  +  1)  where  k  is  the  current 
iteration  number.  If  we  do  this  and  use  his  value  of  b  =  0.1  with  50  bearings, 
we  find  that  RLS  only  converges  when  the  region  of  interest  is  small — so  it 
appears  that  the  values  must  be  chosen  to  reflect  the  typical  distances  between 
the  various  points.  So  a  range  of  6’s  and  A’s  was  tested.  When  b  =  1,  the  final 
estimate  of  emitter  position  oscillates  wildly  for  say  A  <  0.4,  settles  to  within 
the  5cj  uncertainty  ellipse  for  A  =  0.65  — >  0.70,  and  then  falls  quite  short  as  A 
increases:  for  A  =  0.72  the  routine  barely  arrives  at  the  SOu  ellipse. 

As  b  is  increased  over  the  values  10,  10^,  10^  up  to  10®,  the  typical  A  values  for 
which  the  iterations  converge  to  within  the  5a  ellipse  or  better  move  steadily 
upward  through  the  values  0.7,  0.8  etc.,  but  the  estimated  emitter  positions 
still  fall  short  as  A  ^  1,  unless  b  >  10®.  By  the  time  b  >  10^,  then  for  A  >  0.95, 
the  approximation  to  even  2a  or  better  is  very  good — and  it  remains  very  good 
when  the  region  over  which  the  scenario  is  played  out  shrinks  to  the  size  for 
which  Healy’s  values  b  =  0.1,  A  =  1.65k/{k  +  1)  were  usable  above.  Hence  the 
values  chosen  for  the  following  simulations  were  6  =  10®,  A  =  1. 

Further  preliminary  notes  to  emerge  from  running  the  programme  are  summarised  in  the 
following  paragraphs. 


CPLE  approach.  This  gives  quite  accurate  initial  estimates  for  the  case  of  a  stationary 
emitter.  For  example,  if  the  flight  path  containing  ten  data  points  (each  with  2°  bearing 
noise)  subtends  an  angle  of  15°  at  the  emitter,  then  typically  the  emitter  position  is 
calculated  as  fairly  evenly  spread  over  a  2a  Cramer-Rao  ellipse,  with  a  hnal  down-range 
error  of  about  2%,  and  a  final  cross-range  error  of  <1%. 

The  method  also  gives  excellent  estimates  of  initial  position,  velocity  and  acceleration 
for  an  emitter  of  constant  acceleration,  provided  no  more  than  about  0.1°  bearing  noise 
is  present.  Beyond  this  it  loses  accuracy  quickly.  The  more  we  ask  of  our  bearing  data  in 
terms  of  showing  higher  derivatives  of  position,  the  more  sensitive  the  results  are  to  noise. 
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Figure  10:  Scenarios  1  (left)  and  2  (right)  referred  to  on  page  ^ 


Kalman  covariance  matrix.  Although  the  Cartesian  Kalman  covariance  matrix 
represents  the  uncertainty  in  the  state  estimate,  in  practice  its  entries  might  be  huge  even 
when  the  estimates  are  quite  close  to  the  actual  emitter  position,  or  they  might  be  small 
even  though  the  routine  is  not  locating  the  emitter  very  well  at  all.  So  they  seem  not 
to  be  useful,  except  for  the  fact  that  is  used  to  calculate  the  next  state  estimate, 
and  this  can  sometimes  be  very  accurate  even  though  P^j^’s  entries  are  anomalously  large. 
Presumably  the  failure  of  Pk\k  to  adequately  reflect  the  true  situation  is  related  to  the  use 
of  linearisation  in  the  extended  Kalman  filter. 


Cramer— Rao  lower  bound.  This  is  usually,  if  not  always,  much  better  than  the  typ¬ 
ical  accuracies  achieved  by  RLS  and  Kalman,  although  Gauss-Newton  does  tend  to  give 
accuracies  around  the  Cramer-Rao  mark. 


Algorithm  speed  and  number  of  computations.  The  various  algorithms  can  be 
quantified  more  fully,  to  measure  how  long  they  take  and  how  many  computations  they 
each  perform.  Two  simple  scenarios  are  shown  in  Figure  CGI  which  shows  how  these  are 
plotted  within  the  Matlab  gui. 

Scenario  1:  Use  a  stationary  emitter,  and  a  receiver  flying  at  constant  velocity  and  taking 
500  bearings.  Model  the  emitter  as  stationary. 

Scenario  2:  Use  an  emitter  moving  with  constant  speed  but  in  a  slightly  curved  path, 
with  the  receiver  following  a  more  zig-zag  motion.  This  more  complicated  receiver 
motion  must  be  chosen  because  the  receiver  needs  to  outmanoeuvre  the  emitter  in 
order  for  the  geolocation  to  work.  Again  take  500  bearings.  We  try  various  models: 

a:  Emitter  modelled  as  stationary. 

b:  Emitter  modelled  as  accelerating,  with  initial  conditions  supplied  that  are  not 
very  accurate. 

c:  Emitter  modelled  as  accelerating,  with  very  accurate  initial  conditions  supplied. 
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For  these  runs  we  have  used  the  following  parameters  (refer  to  Figure  QH]  for  the  scale): 
initial  position  estimate  =  (600,  700);  initial  velocity  estimate  =  (2,  0);  initial  acceleration 
estimate  =  (0,0);  receiver  x,y  position  errors  taken  to  be  zero  for  all  the  routines,  and: 

CPLE:  Has  no  parameters. 

G— N:  Gain  =  1,  tolerance  =  5  (i.e.  this  is  the  distance  between  subsequent  emitter  posi¬ 
tion  estimates  that,  once  achieved,  causes  the  algorithm  to  stop). 

RLS:  A  =  1,  6  =  10®  as  explained  on  page  [391 

CEKF:  Pq|q  =  1000  x  unit  matrix,  Q  =  unit  matrix. 

MPEKF:  Initial  emitter  distance  error  =  100;  initial  emitter  speed  error  ay  =  5;  initial 
range  ~  500.  This  routine  is  found  to  be  grossly  unstable  if  we  model  the  emitter’s 
acceleration  by  adding  any  sort  of  Q,  so  we  have  set  matrix  Q  equal  to  zero. 

The  performance  of  the  various  routines  is  presented  in  Tables  on  the  next  page.  Each 
main  entry  refers  to  a  1°  bearing  error,  with  the  corresponding  result  for  a  0.1°  bearing 
error  in  parentheses.  What  do  these  results  tell  us? 

CPLE  handles  stationary  emitters  well.  It  handles  accelerating  emitters  well  only  if  they 
are  modelled  correctly,  but  even  then,  is  easily  thrown  out  by  bearing  noise.  It  is 
fairly  fast  to  run. 

Gauss— Newton  takes  the  longest  time  to  process  the  data,  but  this  is  obviously  affected 
by  how  strongly  we  require  it  to  converge,  by  stopping  its  iterations  when  the  final 
result  changes  by  less  than  some  tolerance.  It  locates  stationary  emitters  well.  Like 
the  CPLE,  its  performance  degrades  if  we  model  the  emitter  incorrectly.  Although 
fairly  insensitive  to  our  guess  of  the  emitter’s  initial  position,  it  is  badly  affected  by 
poor  estimates  of  the  initial  velocity  and  acceleration. 

It  should  be  remembered  that  this  G-N  calculation  is  being  done  from  one  position 
of  the  receiver  (the  last),  using  all  of  the  available  data,  and  so  is  an  estimate  made 
at  that  position  only.  On  the  other  hand,  RLS  and  Kalman  give  one  estimate  for 
each  flight  path  point  while  the  receiver  is  moving.  Even  so,  their  estimates  do  not 
always  settle  to  any  obvious  limit  point  at  all.  In  fact  these  estimates  might  wander 
about  signihcantly  before  apparently  beginning  to  settle  in  the  vicinity  of  the  actual 
emitter  position. 

RLS  performs  well  for  a  stationary  emitter  in  the  presence  of  noise.  When  the  emitter 
moves,  good  results  are  only  obtained  if  we  model  it  correctly,  and  then  RLS  is  not 
greatly  affected  by  a  poor  choice  of  initial  conditions.  However  it  can  sometimes 
diverge,  even  for  the  low  bearing  error  case  chosen  of  0.1°. 

Cartesian  EKE  only  performs  well  if  the  emitter  is  modelled  loosely  as  stationary,  al¬ 
though  the  addition  of  a  non-zero  Q  means  that  we  are  actually  modelling  a  certain 
amount  of  unknown  behaviour  in  its  motion.  In  fact,  it  appears  to  be  more  effec¬ 
tive  to  include  a  nonzero  Q  to  model  acceleration,  than  it  is  to  explicitly  model  the 
emitter  as  accelerating  by  altering  F  in  (3.46). 
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Table  1:  Performance  comparison  of  the  different  routines:  Scenario  1 
(page  \4fff-  Main  entry  is  for  1°  bearing  error;  parentheses  refer  to  result 
for  0.1° 


CPLE 

G-N 

RLS 

CEKF 

MPEKF 

Final  cross  range 
error  (%  of  final  range) 
Emitter  placement 
error  (%  of  final  range) 
cpu  time  (s) 
kflops 

0.11  (0.01) 

0.2  (0) 

1  (1) 

23  (23) 

0.04  (0) 

0(0) 

3.7  (3.7) 
81  (81) 

0.3  (0) 

0(0) 

1.5  (1.5) 
50  (50) 

0.12  (0.03)  1.6  (31) 

5  (10)  41  (430) 

1.2  (1.2)  2.5  (2.5) 

59  (59)  433  (433) 

Table  2:  Performance  comparison  of  the  different  routines:  Scenario  2a 
(page  Ifff)-  Main  entry  is  for  1°  bearing  error;  parentheses  refer  to  result 
for  0.1° 

CPLE 

G-N 

RLS 

CEKF 

MPEKF 

Final  cross  range 
error  (%  of  final  range) 

2.3  (2.1) 

2.0  (1.9) 

3(3) 

0.6  (0.1) 

16  (160) 

Emitter  placement 
error  (%  of  final  range) 

55  (56) 

57  (57) 

55  (56) 

8(3) 

50  ±  50  (400) 

cpu  time 

1(1) 

4.5  (4.6) 

1.5  (1.5) 

1.2  (1.2) 

2.5  (2.5) 

kflops 

23  (23) 

99  (99) 

50  (50) 

59  (59) 

433  (433) 

Table  3:  Performance  comparison  of  the  different  routines:  Scenario  2b 
(page  {JW-  (“div.”  =  “diverges” .)  Main  entry  is  for  1°  bearing  error;  paren¬ 
theses  refer  to  result  for  0.1° 


CPLE 

G-N 

RLS 

CEKF 

MPEKF 

Final  cross  range 
error  (%  of  final  range) 

4  (0.05) 

diverges 

(diverges) 

0.3  ±0.2  (0.06) 

div.  (div.) 

div.  (div.) 

Emitter  placement 
error  (%  of  final  range) 

55  (2) 

div.  (div.) 

100  ±30  (1.5) 

div.  (div.) 

div.  (div.) 

cpu  time 

1.2  (1.2) 

div.  (div.) 

2.2  (2.2) 

1.6  (1.6) 

2.5  (2.5) 

kflops 

106  (106) 

div.  (div.) 

480  (480) 

865  (865) 

433  (433) 

Table  4'  Performance  comparison  for  Scenario  2c  (page\40f.  Main  entry 
is  for  1°  bearing  error;  parentheses  refer  to  result  for  0.1° 


CPLE 

G-N 

RLS 

CEKF 

MPEKF 

Final  cross  range 
error  (%  of  final  range) 

4  (0.05) 

0.02  (0.01) 

3  ±  2  (0.02) 

div.  (div.) 

0.7  (0.3) 

Emitter  placement 
error  (%  of  final  range) 

55  (2) 

3.4  (0) 

100  ±  200  (0.3) 

div.  (div.) 

17  (10) 

cpu  time 

1.2  (1.2) 

2.8  (1.6) 

2.3  (2.2) 

1.6  (1.6) 

2.5  (2.5) 

kflops 

106  (106) 

240  (140) 

480  (480) 

865  (865) 

433  (433) 

42 


DSTO-RR-0319 


Table  5:  Performance  comparison  of  the  different  routines 


CPLE 

G-N 

RLS 

GEKF 

MPEKF 

Need  initialising? 

no 

yes 

yes 

yes 

yes 

Correct  convergence  (low  noise)? 

V.  good 

good 

V.  good 

poor 

good 

Robust?  Stationary  emitter: 

yes 

yes 

yes 

no 

no 

Robust?  Moving  emitters: 

yes 

no 

moderate 

no 

no 

Heavily  affected  by  noise? 

yes 

no 

no 

no 

no 

Modified  polar  EKF  This  is  very  sensitive  to  a  good  choice  of  initial  conditions,  as  was 
also  found  by  pT] .  Although  in  principle  capable  of  good  accuracy,  this  routine  needs 
a  lot  of  computation  even  for  the  simple  Scenario  1  on  page  @01  But  certainly  the 
nature  of  the  computations  is  also  important:  for  example,  note  that  in  Scenario  2b 
the  Cartesian  EKF  takes  only  1.6  seconds,  but  uses  many  more  kiloflops  (kilo  floating 
point  operations)  in  Matlab  than  the  other  routines;  those  mostly  take  longer  to  run 
but  use  fewer  kiloflops.  The  main  improvement  of  Modified  Polar  over  Cartesian 
EKF  is  that  MPEKF  is  far  less  sensitive  to  how  we  model  the  emitter’s  motion. 
Good  results  are  obtained  even  if  we  model  a  constant  velocity  emitter  as  having 
constant  acceleration,  whereas  the  CEKF  would  perform  very  badly  in  that  situation. 

In  general,  for  Scenario  1  on  page @01  (stationary  emitter),  CPLE  and  G-N  converge  better 
to  the  actual  position  than  do  either  RLS  or  Kalman.  CPLE  does  not  need  initialising, 
while  G-N  is  extremely  robust  and  will  usually  converge  to  the  correct  position  when 
the  initial  estimate  is  so  bad  that  RLS  and  Kalman  give  useless  results.  Extremely  noisy 
bearings  usually  throw  CPLE,  RLS,  and  Cartesian  Kalman  completely  off  track,  and 
though  they  might  still  converge  to  some  limit  point,  it  will  almost  certainly  be  nowhere 
near  the  actual  emitter  position.  In  the  same  case  G-N  generally  will  still  converge  very 
accurately,  although  it  might  take  longer  than  usual  to  do  this.  Broadly  speaking,  for  a 
one-off  estimate  of  where  the  emitter  is,  G-N  gives  better  results  than  RLS/Kalman,  while 
if  we  are  interested  in  fast  real  time  tracking,  then  RLS/Kalman  are  better  suited  if  they 
are  provided  with  a  good  initial  state  estimate.  This  estimate  can  easily  be  provided  by 
the  CPLE  approach.  These  results  are  summarised  loosely  in  Table  although  it  is  very 
difficult  to  draw  general  conclusions. 


6.1  Reproducing  Some  Known  Results 

In  this  section,  we  run  the  Matlab  programme  to  check  on  some  known  results  and  as  a 
further  comparison  of  the  different  routines.  The  results  we  are  seeking  to  reproduce  are 
samples  taken  from  graphs  in  [34].  The  emitter-receiver  layout  as  used  in  this  reference 
is  shown  in  Figure  TTl  In  the  four  different  scenarios  reproduced  here,  the  receiver  always 
moves  at  105  m/s  westwards  (204  kts)  with  each  flight  lasting  for  350  seconds,  taking  one 
bearing  every  50  seconds.  In  all  cases  the  receiver  position  is  known  with  zero  error. 

Reference  [34]  does  not  specify  the  algorithm  it  uses,  describing  it  only  as  a  simple  least- 
squares  routine.  It  quantifies  this  routine’s  accuracy  for  each  geometry  by  doing  500  runs — 
each  run  giving  one  estimate  of  the  emitter  position — which  are  then  averaged,  producing 


43 


DSTO-RR-0319 


Track  angle 


Midpoint  range 


^ - ! - 

Receiver 

Figure  1 1 :  Typical  geometry  used  in  Section  \6.1\ 


the  mean  deviation  of  the  estimated  emitter  position  from  its  actual  position  (or  its  last 
position  if  it  is  moving).  We  have  done  likewise  using  the  Matlab  programme.  Additionally, 
the  Matlab  programme  produces  a  50%  CEP  that  gives  a  simple  one-parameter  indication 
of  the  spread  of  the  estimates.  In  the  Matlab  programme  each  scenario  uses  only  100  runs, 
since  this  number  already  gives  good  statistics. 

The  four  scenarios  chosen  were  as  follows: 

1.  Emitter  stationary  (at  the  middle  of  the  path  drawn  in  Eigure  TT]);  bearing  error  1°. 

Range  from  midpoint  of  receiver  path  to  emitter  is  50  km. 

2.  Emitter  moves  at  10  m/s  at  a  track  angle  of  30°;  bearing  error  1°. 

Midpoint  range  50  km. 

3.  Emitter  moves  at  20 m/s  at  a  track  angle  of  —60°;  bearing  error  1°. 

Midpoint  range  50  km. 

4.  Emitter  moves  at  20  m/s  at  a  track  angle  of  45°;  bearing  error  2.3°. 

Midpoint  range  200  km. 

Eor  each  scenario  we  first  run  the  CPLE  routine,  always  modelling  the  emitter  as  stationary 
even  when  it  is  moving,  as  per  [34]  •  As  usual,  we  use  CPLE  to  seed  the  other  routines. 
The  only  exception  to  this  is  the  modihed  polar  Kalman  filter,  which  is  fairly  sensitive 
to  initial  conditions,  and  the  forms  it  requires  are  not  quite  given  by  the  CPLE  routine. 
In  particular,  MPEKE  needs  an  initial  speed  error  estimate  which  CPLE  cannot  provide, 
since  we  are  necessarily  modelling  the  emitter  as  stationary.  This  error  has  been  specified 
somewhat  arbitrarily  by  us. 

The  results  for  the  mean  deviation  of  the  estimated  emitter  from  the  actual  emitter 
(or  its  last  position  if  it  is  moving)  are  given  in  Table  [61  The  report’s  values  were  taken 
from  its  graphs,  while  the  values  output  from  the  Matlab  programme  include  a  50%  CEP 
in  parentheses. 

In  general,  the  spread  of  estimates  tends  to  overlap  the  actual  emitter  position.  But 
this  is  by  no  means  always  the  case,  because  of  varying  sensitivities  to  initial  estimates;  in 
particular  the  fact  that  we  are  modelling  the  emitter  as  stationary  means  that  we  are  always 
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Table  6:  Comparison  of  performance  values  of  [Sff  with  the  Matlab  pro¬ 
gramme  described  in  this  report.  Shown  are  results  for  the  mean  deviation  of 
estimated  emitter  position  from  actual  emitter  position  ( or  its  last  position 
if  it  is  moving).  Values  in  parentheses  are  a  50%  CEP.  All  distances  are  in 
kilometres 


Scenario  Report  [34]  Matlab  programme  modelling: 

(page  144|)  (unknown 
algorithm) 


CPLE 

G-N 

RLS 

CEKF 

MPEKF 

1 

1.4 

0.12  (1.1) 

0.14  (1.0) 

0.24  (0.9) 

0.10  (0.8) 

1.4  (4.6) 

2 

1.5 

1.1  (1.4) 

0.9  (1.3) 

0.9  (1.2) 

1.1  (1.1) 

6.1  (5.5) 

3 

12 

11  (11) 

12  (12) 

12  (12) 

13  (13) 

3(6) 

4 

63 

56  (61) 

8  (33) 

47  (36) 

43  (42) 

99  (66) 

estimating  its  velocity  to  be  zero,  and  this  can  badly  degrade  some  routines’  performance. 
Further  simulations  with  the  emitter’s  constant  velocity  being  modelled  cannot  be  done 
in  this  scenario,  because  the  receiver  never  accelerates,  and  so  it  can  never  geolocate  a 
moving  emitter.  Hence  we  have  not  added  any  results  for  the  emitter  modelled  as  moving, 
since  that  requires  a  remodelled  receiver  path,  which  then  complicates  or  destroys  the 
comparison.  But  the  performances  of  the  Matlab  routines  do  generally  agree  with  those 
of  [34],  which  is  what  we  set  out  to  test. 

The  sensitivity  of  MPEKF  to  our  choice  of  initial  conditions  is  worth  noting.  The 
“99  (66)”  values  correspond  to  an  initial  range  of  500  km,  range  error  of  100  km,  and 
speed  error  of  5km/s  (rather  large,  admittedly).  If  we  change  these  values  to  200,  50 
and  0.1  respectively  the  tabulated  values  become  63  (64).  If  we  now  change  the  initial 
conditions  to  200,  10,  0.5,  the  tabulated  values  become  1.3  (46).  These  tabulated  values 
will  also  change  from  one  set  of  100  runs  to  the  next.  Clearly  this  sensitivity  implies 
possible  low  expectations  for  our  results  if  we  start  with  no  prior  knowledge  of  the  initial 
conditions. 


6.2  Checking  an  Old  Suggested  Rule  of  Thumb 

Reference  pM]  writes  down  rules  of  thumb  for  computing  expected  cross-  and  down-range 
errors  for  the  case  of  a  stationary  emitter  and  a  receiver  moving  with  constant  velocity.  It 
includes  no  derivations,  and  its  two  rules  appear  to  be  based  only  on  a  vague  but  incorrect 
geometry.  Certainly  its  expression  for  relative  cross-range  error  is  not  reasonable,  and 
does  not  work  in  practice. 

If  the  cross-  and  down-range  errors  are  labelled  cr^^,  Ojj  respectively,  with  k  bearings 
over  a  flight  path  that  sweeps  out  an  angle  Og  as  seen  from  the  emitter  (as  in  Figure  T2]). 
with  a  bearing  error  of  then  the  rules  of  thumb  are  given  to  within  a  factor  of  about 
unity  as 

^ - ^  ~  /g  .  N 

R  ~  Vk  R~  OsVk'  ^  ’ 


45 


DSTO-RR-0319 


Figure  12:  Typical  scenario  tested  in  SectionlfTM 


Note  that  the  relative  cross-range  error  is  independent  of  swept  path,  something  we  will 
show  to  be  false.  This  is  not  a  simple  mistake  in  }34j.  because  the  relative  cross-range 
error  is  plotted  there  quite  explicitly  as  a  constant. 


Expressions  (16. ip  are  reminiscent  of  Beasley  and  Miles’  result  (14. ip.  which  has 


^  oc  —  ^  oc  ^ 

R  Vk  R  OsVk' 


(6.2) 


However  Beasley  and  Miles’  result  is  calculated  only  for  a  very  specihc  case  of  TDOA  with 
two  receivers,  and  so  cannot  include  any  bearing  error  a^.  Their  cross-range  error  also 
comes  out  as  independent  of  swept  angle  Og- 


The  scenario  used  in  the  present  report  to  test  the  rules  in  (|6.ip  is  shown  in  Figured^ 
A  stationary  emitter  is  located  at  (100,  900),  with  a  receiver  flying  at  constant  velocity, 
starting  at  (100,  100).  The  end  point  of  the  receiver  flight  path  is  placed  variously  over  the 
grid  in  such  a  way  that  the  swept  angle  Og  and  the  angle  ol  are  varied  over  a  convenient 
range.  For  each  of  the  end  points  these  two  numbers  produce,  we  consider  a  range  of 
bearing  errors  and  a  range  of  numbers  of  bearings  taken. 


Each  single  scenario  thus  resulting  is  then  repeated  100  times,  with  new  noise  randomly 
added  in  each  run  for  a  Monte  Carlo  approach.  This  yields  average  measured  values  of 
(Tcr! R  aiid  ^r/ Rj  together  with  the  corresponding  predictions  obtained  from  (6.1). 

As  well  as  the  relative  cross-  and  down-range  errors,  the  circular  error  probable  (CEP) 
is  also  important,  since  such  a  figure  is  necessary  for  deciding  whether  a  missile  can  pass 
close  enough  to  the  emitter  for  its  explosive  radius  to  be  sufficient  to  have  an  effect. 


Comparison  with  Rules  of  Thumb  (“r.o.t.”): 

CPLE  As  the  number  of  bearings  taken  increases,  (7^1  R  increases  from  ~  0.2  x  rule  of 
thumb  to  ~  2  X  r.o.t.,  while  a^^j^jR  increases  from  ~  0.25  r.o.t.  to  5  r.o.t. 
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G-N  UnlR  drops  from  0.15  r.o.t.  to  a  very  small  fraction  if  a  is  small,  but  stays  roughly 
constant  at  0.15  r.o.t.  if  a  is  not  small;  is  always  about  0.20  r.o.t.,  dropping 

slightly  as  the  number  of  bearings  increases. 

RLS  fjfj/ii  drops  from  0.15  r.o.t.  to  <  0.05  r.o.t;  is  always  at  0.30  r.o.t.,  perhaps 

dropping  slightly  with  an  increasing  number  of  bearings. 

CEKF  CTjj/i?  increases  from  0.20  r.o.t.  to  r.o.t;  is  around  0.40  r.o.t.  but  can  fluc¬ 

tuate  wildly. 

MPEKF  (Tfj/i?  drops  from  r.o.t.  to  perhaps  0.1  r.o.t;  (y^^jR  drops  from  0.50  r.o.t.  to 
perhaps  0.05  r.o.t.  if  a  =  0,  but  can  grow  to  2  r.o.t.  if  a  =  60°. 

Clearly  the  trends  are  not  especially  simple  or  linear. 


Least  squares  fitting  of  the  data.  Consider,  say,  a^j^jR  as  a  function  of  fj^,  0g,  a, 
and  the  number  of  bearings  n.  We  wish  to  try  to  ht  some  sort  of  function  of  the  four 
parameters  to  the  data,  in  a  least-squares  sense.  A  cue  for  what  function  to  use  is  taken 
from  the  fact  that  the  data  is  supposedly  at  least  roughly  fitted  by  (j6.1|).  So  perhaps  a 
polynomial  in  a^,9g,a,  and  \/n  is  needed.  This  immediately  presents  a  problem,  because 
even  a  simple  case  can  have  many  terms.  A  “simple”  fitting  function  might  be: 

-h  z^ajd^  -h  . . .  Z2i/Vn  ,  (6.3) 


where  the  highest  order  of  any  parameter  is  just  one,  terms  with  one  parameter  in  nu¬ 
merator  and  denominator  are  included,  and  terms  of  order  —1  are  included  since  we  are 
fitting  a  surface  that  looks  to  have  something  like  a  hyperbolic  shape.  The  calculation  can 
be  done  easily  using  Matlab,  but  a  good  ht  does  not  result:  although  many  data  points 
might  ht  within  a  10%  margin,  others  will  be  out  by  a  factor  of  2.  Of  course  higher  order 
htting  surfaces  can  be  tried,  but  the  number  of  terms  required  to  describe  these  becomes 
very  large.  Even  using  (6.1)  as  a  cue  to  htting  with  just  one  parameter, 


R 


Zl  +  Z2 


(6.4) 


does  not  give  a  good  ht  over  the  whole  surface,  and  experimenting  with  higher  order 
polynomials  in  a^/ {9gy/n)  does  no  better.  So  it’s  not  apparent  that  any  useful  rule  of 
thumb  can  be  arrived  at,  even  with  such  a  simple  geometry  as  tried  here:  a  constant 
velocity  receiver  and  a  stationary  emitter. 


7  Concluding  Remarks 

The  search  for  straightforward  rules  of  thumb  for  expected  geolocation  accuracies  is  a 
difficult  one,  at  least  in  the  sense  of  writing  down  an  expression  that  takes  in  the  param¬ 
eters  of  a  platform’s  motion,  and  produces  the  various  geolocation  errors  for  the  various 
techniques.  The  search  depends  on  which  technique  is  used,  and  also  depends  heavily  on 
receiver  geometry  and  the  kinematics  of  emitter/receiver  motion. 
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In  the  preceding  pages,  we  have  described  tackling  the  problem  by  first  explaining 
common  geolocation  techniques.  Various  details  of  theory  and  implementation  have  been 
given,  so  that  the  algorithms  can  be  run  with  little  extra  research  required  of  the  reader. 

Although  we  have  calculated  some  Cramer-Rao  bounds,  being  the  theoretical  best 
possible  accuracy  obtainable  for  a  given  setup,  we  have  tended  to  take  a  sort  of  “hands-on” 
statistical  view  of  quantifying  the  accuracies.  That  is,  for  this  report  Matlab  software  was 
developed  that  simulates  various  scenarios,  and  this  was  used  in  a  Monte  Carlo  approach 
to  build  up  many  estimates  of  where  an  emitter  is,  based  on  gaussian  noise  added  as 
bearing  error.  Initial  work  was  hampered  by  relatively  long  computing  times,  since  much 
was  being  asked  of  the  software.  Currently,  even  on  a  time  scale  of  one  or  two  years, 
computer  speeds  increase  enough  that  more  and  more  processing  can  readily  be  done  in 
geolocation  algorithms. 

In  the  final  analysis,  perhaps  it  might  be  that  the  best  way  to  predict  the  accuracy 
that  a  geolocation  algorithm  can  achieve,  is  to  model  a  required  scenario  fairly  specihcally. 
In  this  respect,  attempts  to  optimise  a  receiver  configuration  from  the  outset  in  any  purely 
theoretical  way  appear  to  be  very  difficult — especially  without  depending  on  some  sort  of 
visualisation  tool  such  as  described  in  this  report. 
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Appendix  A  Posterior  CRLB 


The  posterior  Cramer-Rao  lower  bound  was  discussed  in  Section  I3.2.11  and  the  ideas 
underlying  it  are  outlined  in  this  appendix. 

When  using  a  recursive  geolocation  routine,  we  need  to  take  into  account  the  prior 
knowledge  available  at  each  step  of  the  recursion.  The  state  x  that  we  seek  now  takes  on 
a  new  value  as  each  measurement  arrives.  Hence,  index  it  with  the  measurement  number: 
the  measurement  set  . . . ,  Zp,}  estimates  the  latest  value  of  the  state,  Prior 

information  is  given  by  the  quantity  Zg,  which  will  be  a  measurement,  but  is  not  counted 
as  part  of  the  data  set. 

Instead  of  dealing  with  the  conditional  density  a; as  in  Section  13. 2.11  the  state  Xp. 
is  now  treated  as  a  random  variable  in  its  own  right,  and  so  the  combined  density  p{Zp,,Xp,) 
is  used  to  compute  the  bayesian  Fisher  information  at  each  time  k: 

Jb(I=)  s  (vLy(Zt,%))}  ,  (Al) 

where  Define  L{xp.)  via 

L{Zk,x,^)  =  -lnp(Zfc,a;fc) 

=  -lnp(a;fe)  -lnp(Zfc|a;fc) 

=  L(a;^)  +  L(Z^|a;^) .  (A2) 


Equations  (]All  and  (]A2^  combine  to  give 


^ - V - ^ 

=  Jp(k),  the  prior  information 


(vLy(z»k.))}. 

^ - V - ^ 

=  data  information 


(A3) 


where  we  have  replaced  by  E^,^  in  the  first  term  on  the  right  in  (1A3) ,  since  its 

argument  does  not  depend  on  the  measurements  Zp,.  Also,  as  is  typically  done  to  simplify 
calculations  in  bayesian  theory,  the  second  term  of  (A3) .  E^,  usually  replaced 

by  Jj\i^{k).  This  corresponds  to  using  an  unchanging  set  of  initial  conditions  in  a  set  of 
Monte  Carlo  simulations. 

Suppose  that  an  initial  measurement  Zg  (i.e.  prior  information)  has  indicated  a  target 
state  Xq,  and  assume  Xq  is  normally  distributed:  p{xq)  ~  M{xq,  Pq),  by  which  is  meant 

„  I  exp  {xq  -  Xof  Pg-i  {xq  -  Sg)  .  (A4) 

V  I  ^Ti'Po  I  ^ 

Then  if  Up  and  Vp  are  set  to  zero  in  (3.46)  (as  done  in  bearings-only  tracking),  it  follows 
that  Xp  is  normally  distributed:  p{xp)  ~  M{xp,  Pp).  Setting  Fp  =  F,  a  constant  matrix, 
the  relevant  expressions  now  involve  the  power  of  P: 

Xp  =  F^x,,  P,  =  P"Po(P^)^  (A5) 


In  that  case,  the  prior  information  turns  out  to  be 

Jpik)=Pp\ 


(A6) 
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Thus,  the  calculation  of  the  posterior  CRLB  in  (]A3f)  reduces  to  a  calculation  of  Jj^^,  the 
non-bayesian  Fisher  information.  The  calculation  of  this  matrix  is  somewhat  laborious; 
the  remainder  of  this  appendix  outlines  the  way  it  is  done,  for  an  emitter  modelled  as 
having  constant  velocity. 

For  this  constant  velocity  case,  the  emitter  state  is  defined  in  terms  of  its  position 
i^k-’Vk)  velocity  {xk,yk)  at  time  k: 

^k=[xk  Vk  ijk]^  ■  (A7) 

Beginning  with  (j3.3[  13. 6j) ,  we  need  the  conditional  probability 


P{Zk\Xk) 


y/2'KR 


exp 


-1 


k 


h{Xk,i)f  , 


(A8) 


which  directly  gives  L{Z^\x^).  The  quantity  h{Xf,,i)  is  the  same  as  used  in  (j3.52[).  being 
a  sort  of  backwards  prediction  of  the  exact  bearing  that  should  be  seen  at  time  i  ^  k, 
given  the  current  time  k  and  current  state  x^,.  With  equal  time  steps  At,  we  can  write 
the  emitter’s  position  at  some  intermediate  time  i  as: 


ixi,yi)  =  {xk,yk)  + 


(A9) 


so  that  (omitting  the  quadrant-dependent  constant  in  the  following  inverse  tan  functions, 
which  will  differentiate  to  zero  anyway — see  the  comment  on  page  |T9|) 


h{x^,i) 


tan-i  =  tan-i  ^ 

Xi  -  Axj 

-1  +  A,, 


tan 


(for  all  emitter  models) 

(constant  velocity  model  only),  (AlO) 


where  is  the  receiver  position  at  time  i.  The  Fisher  information  matrix  has  the 

raw  form  of 

(vi^L(Z,lx,)^  }  .  (All) 

For  convenience,  if  we  rewrite  the  state  as 


r  1  2 

[Xk  x^ 


(A12) 


then,  with  further  effort,  the  Fisher  information  matrix  can  be  shown  to  have  as  its  (m,  n)*^ 
element; 


Jnsik) 


mn 


1  s;^dh{Xf^,i)  dh{x^,i) 
R  ^  dx^  dx'^ 


(A13) 


Actually  this  last  expression  holds  for  emitter  models  with  arbitrary  motion,  not  just 
constant  velocity;  what  does  change  for  different  emitter  models  is  the  form  of  h{x^,i) 
in  the  second  line  of  (AlO).  where  the  numerator  and  denominator  would  contain  extra 
Taylor  series  terms  involving  higher  derivatives  of  position. 
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For  the  constant  velocity  case,  further  manipulation  leads  to 


1 


E 


1 


Ax? 


(symmetric) 


Ay?(i  -  /c)At 
-A?/-Axj(i  -  k)M 


-Ay-Axj(i  -  k)At  ' 
Ax?(i  —  k)At 
—Ay-Ax^{i  —  k)‘^At^ 
Ax^(z  —  k)‘^At‘^_ 

(A14) 


The  bayesian  Fisher  information  J^{k)  then  follows  from  (A3[  A6i  A14). 
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Appendix  B  Total  Least  Squares 


Here  we  supply  the  relevant  details  referred  to  in  Section  13.6.11 


When  the  matrix  H  of  (|3.14[)  itself  contains  noise,  the  usual  pseudo  inverse  is  not 
guaranteed  to  give  the  most  accurate  estimate  of  the  state  x.  In  this  case  the  method  of 
Total  Least  Squares  can  sometimes  give  better  results.  To  construct  the  method,  begin  by 
rewriting  (’13.141  as 


[H  -z] 


X 

1 


residual. 


(Bl) 


Suppose  we  multiply  both  sides  by  [H 


t 


[H  -z]‘  [H  -z] 


=  \H  —z]*  residual. 


(B2) 


In  a  noiseless  world  the  right  hand  side  would  be  zero,  in  which  case  A  =  \JL  —2]  *  [if  — z] 
(a  square  matrix)  would  be  singular.  In  practice  we  will  perturb  it  to  produce  a  new 
singular  version,  in  the  least  perturbative  way,  by  zeroing  its  smallest  eigenvalue.  This 
can  be  done  through  first  finding  three  matrices  U,  S,  V  that  make  up  the  singular  value 
decomposition  of  [if  — z] : 

[H  -z]  =  USV^ 

Matlab:  [U,  S,  V]  =  svd([H,  -z]  ,  0).  (B3) 


The  matrices  U,  S  and  V  have  useful  properties  that  need  not  concern  us  here,  except 
that  S  is  diagonal,  having  as  its  entries  the  square  roots  of  the  eigenvalues  of  A.  If  we 
arrange  these  in  descending  order  down  the  diagonal]^  then  we  can  perturb  A  to  be  singular 
by  zeroing  the  last  eigenvalue  of  S,  to  produce  a  new  matrix  S.  In  that  case,  the  new 
matrix  USV^  defines  new  H  and  z  matrices: 


[H  -z]  =  USV* . 


(B4) 


Solving  the  new  system 


[H  -z] 


X 


=  0  ,  i.e.  z  =  Hx  , 


(B5) 


produces  a  possibly  better  estimate  of  the  state  x.  In  practice,  the  calculation  is  further 
simplified  because  the  equation  which  is  to  be  solved. 


USV* 


=  0. 


(B6) 


can  be  rewritten  in  terms  of  the  columns  of  the  n  x  n  matrix  V.  If  the  state  x  has  n 
entries,  then 


for  each  row  i. 


n— 1 


^  [  t/jj  Sjj  Vcolumn  j  —  0  ) 

i=i 


(B7) 


and  V  would  have  to  be  altered  to  suit,  but  this  preferred  arrangement  is  automatically  produced 
by  Matlab  anyway.  We  need  do  nothing  extra. 
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so  that 


X 


must  be  orthogonal  to  the  first  n  —  1  columns  of  V.  But  a  special  property 


of  V  is  that  it  is  orthogonal:  its  columns  are  all  mutually  orthonormal.  Hence 


must 


be  proportional  to  the  last  column  of  V.  In  that  case,  the  total  least  squares  estimate  of  x 
is  iust 

■  Vln 


X  = 


1 

1^ 


K- 


n— l,n 


(B8) 


This  is  all  easily  programmed  in  Matlab,  especially  since  the  matrices  U,  S  are  not  used. 
In  summary,  when  solving  for  x  with  a  noisy  H: 


Hx  =  z  +  noise, 


(B9) 


the  two  approaches  of  least  squares  and  total  least  squares  can  be  programmed  in  the 
following  way: 


Least  squares: 


X  =  H*z 

Matlab:  x  =  pinv(H)  *  z 


(BIO) 


Total  least  squares: 


USV*  =  \H  —z]  (s.v.d.);  x  = 

^nn 


Vln 


Vn—l,n 


Matlab:  [notNeeded,  notNeeded,  V]  =  svd([H,  -z]  ,  0) ; 

X  =  V(l: end-1,  end) /V (end,  end) 


(Bll) 
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Appendix  C  Running  the  Matlab  Geolocation 

Simulator 

The  Matlab  programme  referred  to  in  this  report  is  used  to  investigate  the  tracking  al¬ 
gorithms  of  Sections  I3.7h3.10l  It  allows  us  to  specify  fairly  general  flight  paths  for  one 
receiver  and  one  emitter.  It  then  generates  noisy  bearings  at  regular  intervals  along  the 
flight  paths,  and  uses  our  choice  of  either  the  CPLE,  Gauss-Newton,  RLS,  or  either  type  of 
Kalman  filter  to  estimate  the  emitter’s  flight  path.  It  also  plots  Cramer-Rao  uncertainty 
ellipses  as  described  in  Section  15.1  and  [4] ,  that  serve  as  a  guide  to  how  accurately  we  can 
ever  hope  to  locate  the  emitter. 

Here  we  describe  two  possible  ways  of  running  the  Matlab  gui. 

1.  To  run  the  gui  in  its  simplest  form,  type  geolocationGui  at  the  Matlab  command 
prompt,  which  will  run  the  geolocationGui  .m  file.  This  brings  up  the  gui.  Param¬ 
eters  can  then  be  entered  into  various  boxes.  Pressing  the  appropriate  button  draws 
and  processes  the  flight  paths. 

2.  If  a  better  feel  for  the  statistics  of  the  final  estimation  is  needed,  a  Monte  Carlo 
approach  can  be  used.  This  runs  the  relevant  routine  many  times  for  various  emitter- 
receiver  geometries.  Instead  of  typing  geolocationGui  from  the  command  line,  type 
manyRuns  instead,  which  runs  the  code  in  manyRuns  .m.  The  gui  will  start  as  usual,  at 
which  point  the  appropriate  parameters  must  be  entered  into  it;  this  is  just  equivalent 
to  running  geolocationGui  .m  for  each  set  of  parameters  and  entering  these  by  hand 
into  the  gui.  Currently  many  Runs,  m  is  only  set  up  to  handle  a  stationary  emitter 
and  constant  velocity  receiver. 


Running  geolocationGui  .m  (the  Gui) 


After  typing  geolocationGui,  the  window  in  Figure  Cl  appears.  The  plot  area  for 
flight  paths  is  a  flat  area  with  north  upwards  and  east  to  the  right.  It  starts  off  with  a 
nominal  axes  scaling,  although  the  scaling  can  be  changed  at  any  time  by  clicking  the 
Change  axes  limits  button. 


Once  appropriate  axes  limits  have  been  chosen,  the  paths  of  both  emitter  and  receiver 
are  plotted  on  this  map.  They  can  be  plotted  in  either  order,  and  are  both  recorded  in 
the  same  way.  For  example,  we  enter  emitter  waypoints  by  first  clicking  on  the  button 
labelled  Start  recording  emitter  waypoints  ,  then  clicking  the  left  mouse  button  to  place 
each  waypoint  on  the  map.  We  finish  by  clicking  the  same  gui  button,  now  relabelled 
Stop  recording  emitter  waypoints  . 


The  actual  flight  speeds  and  bearing  measurement  times  involved  are  changed  by  en¬ 
tering  numbers  in  the  two  input  boxes:  Time  between  receiver  waypoints  and  Time 
between  bearing  measurements.  The  time  between  emitter  waypoints  cannot  be  set, 
since  it  is  determined  by  the  number  of  waypoints  for  emitter  and  receiver  we  input, 
together  with  the  time  between  receiver  waypoints. 
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Error  specifier 


Routine  chooser 


Track  recorder 


Monte  Carlo 
specifier 

J 


Initial  conditions 
specifier 


Results  panel 


Figure  Cl:  Startup  window  for  the  geolocation  gui 


As  an  alternative  to  specifying  the  flight  paths  t 
approach.  First  click  on  the 


lis  way,  we  can  use  an  initial  condition 
Specify  initial  conds  button.  This  brings  up  a  new  set  of 
buttons  and  input  boxes.  Before  doing  anything  else,  ensure  that  the  time  between  bearing 
measurements  has  been  entered.  Then,  for  each  of  the  receiver  and  emitter  (in  either 
order),  enter  the  flight  time  (which  should  be  the  same  for  both),  and  the  initial  position, 
velocity,  etc.  Clicking  on  Draw  emitter  path  or  Draw  receiver  path  then  plots  a  path 
based  on  these  initial  conditions. 


Next  specify  how  many  runs  we  wish  to  make.  Specifying  one  run  makes  the  receiver 
and  emitter  fly  their  paths  just  once;  the  receiver  takes  noisy  bearing  measurements  at 
intervals  given  by  the  Time  between  bearing  measurements  input  box,  and  an  answer 
for  the  estimated  position  or  track  of  the  emitter  is  computed.  If  a  number  greater  than 
one  is  chosen,  this  process  is  repeated  many  times,  for  a  new  set  of  noise  values  each  time. 
Each  estimated  emitter  track  is  overlaid  on  the  screen  at  the  end  of  the  computation — 
or  whenever  Matlab  decides  to  draw  its  graphics  buffer  to  the  screen,  which  currently  is 
always  at  the  end  of  the  whole  set  of  computations. 


Checking  the  Retain  estimated  emitter  paths  box  means  that  no  matter  how  many 
runs  are  made,  successive  paths  will  remain  on  the  screen.  To  clear  these  runs,  click  the 
Clear  estimates  button.  This  prevents  Matlab  from  doing  a  slow  screen  refresh  with 
these  paths  in  place  as  part  of  its  screen  management  before  it  starts  the  next  run. 
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Now  enter  the  bearing  error  into  the  Bearing  error  (degrees)  box.  (Note  that  all  of 
the  input  boxes  are  independent  of  each  other,  so  the  various  inputs  can  be  entered  in  any 
order.)  This  is  a  value  that  will  become  the  standard  deviation  of  a  normal  distribution 
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centred  about  0°,  from  which  Matlab  will  draw  random  values  to  add  to  the  exact,  true, 
bearing.  This  noisy  bearing  is  what  Matlab  then  actually  uses  to  perform  its  geolocation 
computation. 

We  can  also  enter  receiver  position  errors  into  the  Receiver  x  error  and  Receiver  y 
error  boxes.  These  numbers,  analogously  to  the  bearing  error,  are  the  standard  deviations 
of  normal  distributions  about  zero  from  which  will  be  randomly  drawn  numbers  to  add 
to  the  receiver  x  and  y  positions.  They  reflect  our  uncertainty  of  where  the  receiver  is, 
due  to  INS  and  GPS  errors.  These  numbers  are  ordinarily  quite  small,  if  they  are  present 
at  all. 

A  geolocation  method  must  be  selected  from  the  following: 


CPLE.  To  run  this  routine  first  click  the  CPLE  button,  which  causes  a  new  set  of 
buttons  to  be  drawn.  There  is  just  one  button  to  press  to  run  the  routine:  (another) 
CPLE  ,  since  this  routine  requires  no  seed  estimate.  Pressing  this  button  causes 
Matlab  to  do  the  following  with  the  flight  paths.  Eirst  it  interpolates  the  waypoints 
so  as  to  make  the  same  number  of  points  in  each  track,  spaced  in  time  by  the  bearing 
measurement  interval.  That  way  the  new  waypoints  can  be  paired  off,  so  that  the 
emitter  position  corresponding  to  each  receiver  measurement  is  unambiguous.  Then 
a  cubic  spline  curve  is  computed  through  each  set  of  points  to  make  two  smooth 
paths. 


Next,  Matlab  builds  a  list  of  noisy  bearings,  one  for  each  receiver-emitter  waypoint 
pair.  Once  this  is  done,  further  processing  of  the  noisy  bearings  depends  on  the 
routine. 


If  we  have  modelled  the  emitter’s  motion  as  constant  position,  then  the  estimate 
produced  by  the  CPLE  routine  will  be  a  single  point.  Likewise  if  we  have  chosen 
constant  velocity  as  the  model,  the  routine  works  out  an  initial  position  and  initial 
velocity,  and  plots  the  constant  velocity  path  resulting  from  these,  which  is  neces¬ 
sarily  a  straight  line.  (This  fact  should  be  noted  because  the  other  routines  might 
do  this  differently;  e.g.  the  Kalman  routine  tracks  the  emitter,  and  will  seldom — 
essentially  never — produce  a  straight  line  for  the  same  constant  velocity  emitter.) 
Similarly,  a  set  of  initial  conditions  is  produced  and  used  to  plot  the  emitter  motion, 
when  we  model  that  motion  as  having  higher  derivatives. 


After  running  the  CPLE  routine,  its  results  can,  if  we  choose,  be  used  as  initial 
conditions  of  any  of  the  other  routines,  which  all  require  a  seed.  Do  this  by  clicking 
the  Use  result  to  set  new  ICs  button. 


Gauss— Newton.  Choose  this  with  the  Gauss-Newton  button.  Enter  a  gain  and  tol¬ 
erance  in  the  box  that  appears.  The  gain  is  a  factor  around  unity  that  multiplies 
the  stepping  size  calculated  by  Gauss-Newton,  in  the  hope  of  helping  the  iterations 
converge  faster.  The  tolerance  is  the  distance  on  the  plot  area  such  that  once  the 
change  in  (hnal)  position  from  one  iteration  to  the  next  is  less  than  this  amount,  the 
routine  will  terminate. 


Next,  estimate  the  initial  conditions.  They  might  have  been  already  estimated  by  a 
preliminary  run  of  the  CPLE;  but  if  not,  choose  an  emitter  motion  model  from  the 
list  next  to  Emitter  is  modelled  with  constant:.  It  is  possible  to  “cheat”  by 
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hitting  the  Guess  ICs  button.  This  calculates  the  actual  initial  conditions  for  the 


emitter  path  that  was  entered,  and  is  a  way  of  entering  something  at  least  half  way 
reasonable,  as  a  way  of  experimenting  with  different  scenarios  for  which  guessing 
reasonable  initial  conditions  is  difficult. 


Now  click  the  newly  appeared  Gauss-Newton 


button  to  run  the  routine.  Its  it¬ 
erations  stop  once  the  tolerance  has  been  reached.  Iterations  also  stop  once  25  of 
them  have  occurred,  since  this  indicates  that  the  routine  is  probably  not  converging 
and  should  be  abandoned.  (This  number  25  is  just  the  value  of  the  maxiterations 
variable  in  Gauss — Newton. m  and  can  be  changed  there  if  required.) 

As  with  the  GPLE  routine,  Gauss-Newton  calculates  a  set  of  initial  conditions  for 
the  emitter.  Thus  if  the  emitter  is  modelled  with  constant  velocity,  an  absolutely 
straight  line  will  necessarily  be  drawn. 

Recursive  Least  Squares.  Glicking  Recursive  least  squares  prints  the  relevant  but¬ 
tons  to  the  gui.  Initial  conditions  and  a  model  of  the  emitter’s  motion  must  be 
specified. 

There  are  now  two  parameters  to  specify  to  run  RLS:  Lambda  and  Q,  as  defined  in 
Section  f3.9i  Typical  values  might  be  Lambda  =  1  and  Q  =  le8  (i.e.  10®).  Pressing 


the  new  Recursive  least  squares 


button  runs  the  routine. 


The  type  of  track  to  expect  of  the  estimated  emitter  plot  is  different  for  this  and 
the  next  techniques.  For  example,  if  we  choose  the  emitter  as  being  modelled  with 
constant  position,  then  unlike  the  GPLE  and  G-N  case,  the  resulting  plot  will  almost 
certainly  not  be  just  one  point.  The  reason  is  that  RLS  and  the  Kalman  filters  are 
not  batch  processors,  and  so  do  not  produce  just  one  state  estimate  from  their 
processing.  Rather,  they  plot  a  point  whose  location  is  then  updated  as  each  new 
bearing  arrives.  Similarly  if  we  choose  to  model  the  velocity  as  constant:  a  straight 
line  is  not  expected  to  be  output,  since  as  before,  the  position  of  the  emitter  is  simply 
being  estimated  and  updated  for  each  bearing. 


Cartesian  Extended  Kalman.  Glick  on  the  Gart.  Kalman  button  to  run  this.  Besides 
our  supplying  initial  conditions  as  before,  the  routine  also  requires  the  parameters 
labelled  P_0 1  0  and  Q  as  defined  in  Section  [3. lOi  Typical  values  might  be  1000  and  1. 
Then  click  the  Gartesian  Kalman  button. 


Modified  Polar  Extended  Kalman.  For  this,  press  the  MP  Kalman 
ter  values  for  sigma_r,  sigma_v  and  Approx,  range 


button  and  en- 


These  are,  respectively,  the 
errors  in  the  initial  position  and  speed,  and  an  estimate  of  the  range  (some  interme¬ 
diate  value  can  be  given  if  the  range  changes  greatly  for  the  supplied  paths).  Then 
click  on  the  Modified  Polar  Kalman  button. 


Other  Buttons  and  Ontput  Fields 

Display  varions...  There  is  sometimes  a  need  to  produce  various  numbers  indicating 


the  performance  of  an  algorithm.  Glicking  on  this  button  writes  these  numbers  to  the 
Matlab  command  line.  It  will  produce  an  error  message  when  clicked  if  some  of  the 
information  coded  to  be  output  doesn’t  exist. 


60 


DSTO-RR-0319 


Axes  ^  eps  Clicking  this  produces  an  encapsulated  postscript  file  of  the  current  plot. 
The  allocated  filename  is  time-tagged  to  prevent  accidental  overwriting. 


Reset  This  clears  all  tracks  from  the  gui  and  sets  various  other  quantities  to  sensible 
values.  Use  this  when  a  new  set  of  flight  paths  is  required  to  be  drawn.  If  we  just  add 
more  waypoints  to  a  plot  without  hitting 


Reset  first,  they  will  be  added  to  the  existing 


flight  paths. 


Emitter  position  estimation  This  button  implements  the  work  of  Kim  Brown  as 
outlined  in  Section  f3.2i  Clicking  on  it  brings  up  a  separate  window  for  parameter  input. 
The  scenario  it  is  considering,  described  comprehensively  in  [4],  is  that  of  a  plane  flying  a 
great  circle  over  Earth,  taking  bearings  of  a  stationary  emitter.  The  numbers  to  be  input 
relate  to  Figure  IC2i 


Broadside  angle  (degrees)  As  in  Figure  IC2i 

Emitter  range  (km)  The  units  of  km  are  purely  arbitrary. 

Bearing  error  (degrees) 

Number  of  points  This  is  the  number  of  points  along  the  receiver  flight  path,  at  each  of 
which  the  receiver  takes  a  bearing  measurement. 

Receiver  speed  (km/h) 

Time  (s)  between  bearing  measurements  The  unit  of  seconds  is  arbitrary.  The  length 
of  the  receiver  path  flown  depends  on  this  time  and  the  number  of  points  specified 
along  the  path. 


Clicking  the  Emitter  position  estimation 
in  such  a  way  as  to  draw  an  uncertainty  e’ 


button  calculates  the  Cramer-Rao  lower  bound 
llipse,  which  is  really  meant  to  be  superimposed 
on  the  emitter  position.  This  is  a  measure  of  the  most  optimistic  error  we  can  hope  for  in 
the  specified  situation. 


This  emitter  position  estimation  theory  has  also  been  implemented  in  the  main  geolo¬ 
cation  routines,  provided  we  are  modelling  the  emitter  as  stationary  there.  If  we  are,  then 
the  same  theory  is  used  to  draw  a  Cramer-Rao  error  ellipse,  this  time  superimposed  on 
the  actual  position  of  the  emitter.  This  is  useful,  because  it  allows  us  to  compare  the  ge- 
olocated  position  with  the  Cramer-Rao  error  ellipse.  However,  note  that  the  Cramer-Rao 
theory  applies  to  unbiased  techniques,  and  the  techniques  used  in  the  gui  are  not  neces¬ 
sarily  unbiased — as  is  demonstrated  whenever  a  spread  of  Monte  Carlo  emitter  estimates 
does  not  cluster  about  the  true  emitter  position. 


Plotting  the  ellipse  is  governed  by  an  edit  field  in  the  main  gui:  What-sigma  ellipses 
to  plot.  If  set  to  e.g.  2,  then  the  ellipse  resulting  is  really  a  contour  of  the  probability 
distribution  at  two  standard  deviations.  For  accurate  geolocation  scenarios,  the  la  ellipse 
is  usually  too  small  to  see,  so  raising  the  value  of  a  is  equivalent  to  magnifying  the  ellipse 
to  make  it  visible.  The  Draw  ellipses  button  must  be  checked  to  switch  this  facility  on, 
since  plotting  the  ellipses  can  be  computationally  expensive. 
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North 


Figure  C2:  Emitter  position  estimation  scenario 


Display  tracks  slowly  Checking  this  button  draws  the  tracks  point  by  point,  as  a 
visual  aid  only.  It’s  only  enabled  for  the  recursive  routines,  to  emphasise  the  fact  that 
these  routines  operate  by  updating  the  emitter  estimate  as  each  bearing-capture  point 
of  the  flight  path  is  reached.  Depending  on  which  routine  has  been  chosen,  only  certain 


combinations  of  this  button  together  with  the 
a  coding  choice,  due  to  the  practicalities  of  ca 
for  the  recursive  routines. 


Draw  ellipses  button  are  enabled.  This  is 
culating  the  posterior  Cramer-Rao  bound 


How  many  runs?  If  we  want  to  do  a  Monte  Carlo  run,  enter  a  number  greater  than  1 
here.  In  principle,  after  each  run  the  resulting  estimate  of  the  emitter  path  will  be  drawn 
to  the  screen,  as  well  as  the  current  run  number  being  printed  on  the  Matlab  command 
line;  but  in  practice,  Matlab  only  draws  its  graphical  buffer  at  the  end  of  the  entire  Monte 
Carlo  set.  Next  to  this  box,  two  values  for  the  circular  error  probable  will  be  printed, 
for  50%  and  95%.  That  is,  the  50%  number  is  the  distance  out  from  the  emitter’s  final 
position  that  we  must  go  in  order  to  encompass  50%  of  the  estimates  to  that  position. 
Similarly  for  the  95%  number. 


Here’s  how  the  routine  performed  After  every  run,  several  performance  indicators 
are  printed  into  this  rectangle  in  the  gui.  These  are: 

Angle  subtended  by  the  emitter  path  in  receiver  frame  (degrees)  This  only  has 
real  use  if  the  receiver  is  not  following  an  overly  convoluted  path.  The  label  is  worded 
in  this  way  to  take  into  account  the  fact  that  the  emitter  might  be  moving.  In  the 
simplest  case  of  a  stationary  emitter,  this  angle  is  identical  to  the  angle  subtended 
by  the  receiver  path  at  the  emitter. 
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Final  rel.  cross-range  error  Refer  to  Figure  [2]  for  the  definition  of  cross-range  error. 
The  relative  cross-range  error  is  the  cross-range  error  divided  by  the  range.  This 
quantity  is  really  the  answer  to  the  more  specihc  question:  at  the  hnal  position  of 
the  receiver,  what  is  the  relative  cross-range  error  in  the  final  estimated  position  of 
the  emitter? 

Final  rel.  down-range  error  This  is  analogous  to  the  relative  cross-range  error. 

Elapsed  cpu  time  for  this  simulation  (s)  The  elapsed  cpu  time  in  seconds  is  a  good 
parameter  to  compare  the  speed  of  the  different  routines. 

Centroid  displacement  from  true  position  Whether  we  are  running  in  Monte  Carlo 
mode  (in  which  case  a  set  of  emitter  paths  will  be  plotted),  or  not  (in  which  case  just 
one  emitter  path  will  be  plotted),  there  will  be  one  or  more  estimates  for  the  final 
emitter  position.  This  set  has  a  centroid.  The  displacement  of  this  centroid  from 
the  true  emitter  position  is  a  good  indicator  of  how  well  the  geolocation  routine  is 
performing. 

In  the  case  of  running  Gauss-Newton,  the  centroid  placement  might  be  cause  for 
confusion  in  that  one  might  ask:  why  is  the  emitter  placement  error  (output  after 
the  calculation)  sometimes  greater  than  the  tolerance  that  we  have  specihed?  This 
is  because:  in  specifying  the  tolerance,  we  are  telling  the  algorithm  to  stop  when  the 
distance  between  successive  estimates  to  the  emitter  position  becomes  smaller  than 
the  tolerance,  so  that  it  doesn’t  keep  iterating  unnecessarily  for  no  real  improvement 
in  its  estimate.  On  the  other  hand,  the  emitter  placement  error  is  the  distance 
between  the  last  point  and  the  actual  emitter  position.  It  might  be  that  G-N  has 
done  a  bad  job,  and  converged  on  a  point  far  from  the  actual  emitter  position.  In 
that  case  the  emitter  placement  error  will  be  large,  even  though  the  iteration  steps 
have  reduced  to  the  step  size  required  to  stop  the  algorithm. 

Running  manyRuns.m  (Batch  Mode) 

The  file  manyRuns.m  contains  code  that  runs  a  batch  of  geolocation  scenarios.  The  basic 
layout  is  covered  in  Section  [6.21  and  shown  in  Figure  [121  on  page  H6l  The  emitter  is  at  rest 
on  the  y-axis,  while  the  emitter  has  a  constant  velocity  at  some  angle  a  to  the  x-axis.  It 
sweeps  out  a  path  of  angle  Og  as  seen  by  the  emitter,  in  a  time  At,  and  it  makes  bearing 
observations  at  the  rate  of  one  per  second.  The  code  is  run  by  typing  manyRuns  at  the 
Matlab  prompt,  at  which  we  then  supply  values  for  the  various  parameters  used  to  control 
the  runs,  as  arrays.  These  are  the  bearing  error,  Og,  a,  and  At.  So  to  enter  an  array  of 
values  0,  5, 10, 15, . . . ,  80,  we  just  enter  0:5: 80,  while  to  enter  an  array  with  a  changing 
increment,  we  enter  the  array  explicitly:  [2  3  5  7  11] 

The  number  of  Monte  Garlo  runs  and  the  routine  used  must  also  be  specified  at  the 
prompt.  All  runs,  apart  from  when  using  the  GPLE  routine,  are  automatically  seeded  by 
hrst  running  GPLE  (which  needs  no  initial  conditions).  The  chosen  routine  is  then  run  as 
many  times  as  entered  for  the  Monte  Garlo  number. 

manyRuns.m  produces  two  files.  The  main  one  is  a  .mat  file  whose  name  we  supply 
when  prompted  by  manyRuns.m.  This  name  is  time-tagged  to  prevent  any  accidental 
write-over:  if  we  specify  the  name  fred  at  time  2:02:41  pm,  the  actual  name  will  be 


63 


DSTO-RR-0319 


fred_23_Mar_2005_14_02_41  .mat 

The  second  file  produced  is  a  “diary  file”  fred_23_Mar_2005_14_02_41  .diary  containing 
a  history  of  the  runs,  in  case  this  is  needed  for  checking  in  the  event  of  a  computer  crash; 
this  can  be  deleted  when  manyRuns.m  finishes. 

The  .mat  file  is  then  read  by  PlotMatrixCurves .m,  which  reads  its  data  and  plots 
various  graphs.  (PlotMatrixCombined .  m  also  does  a  similar  job,  but  is  a  little  less  evolved. 

It  puts  all  of  its  plots  onto  one  page,  in  the  event  that  they  are  to  be  included  this  way  in 
a  report.)  PlotMatrixCurves  .m  asks  for  several  inputs: 

Name  of  input  file.  That  is,  the  above  .mat  file  containing  the  plot  data. 

What  to  plot  on  y-axis.  Refer  to  Section  16.21  Choices  are  relative  down-  and  cross¬ 
range  errors,  and  circular  error  probable.  Choosing  e.g.  relative  down-range  error, 
the  relevant  choice  is  presented  as  “Relative  down-range  error  and  sigma_R/R”. 
What  this  indicates  is  that  the  relative  down-range  error  (the  statistical  error  output 
from  geolocationGui  .m  run  one  or  more  times  through  manyRuns.m)  will  be  plotted 
on  the  same  graph  as  the  “rule  of  thumb  error”  from  (j6.1). 

Bearing  error  values.  Here  we  are  presented  with  a  range  of  bearings  (in  degrees)  that 
have  been  used  to  create  the  data  now  held  in  the  .mat  file.  For  example,  a  line 
saying  0:0. 1:0. 5  means  that  all  the  values  0,  0.1,  0.2,  ...,  0.5  have  been  used  as 
bearing  errors.  Either  enter  one  of  these,  or  enter  nothing  at  all  (just  hit  return). 
Entering  a  specific  number  means  we  wish  to  plot  a  slice  through  the  matrix  of  data 
at  this  value  of  bearing  error.  If  nothing  is  entered,  all  of  the  bearing  errors  will 
become  the  x-axis.  Note  that  if  only  ever  one  bearing  error  was  already  stipulated 
in  manyRuns.m,  then  we  will  not  be  asked  to  specify  that  bearing  again. 

theta_s  values.  Similarly  we  enter  a  value  for  the  swept  angle  9s  in  degrees.  Note  that 
if  we  enter  nothing  here  as  well  as  for  the  bearing  error,  PlotMatrixCurves  .m  will 
crash  as  it  tries  to  put  two  sets  of  values  along  the  x-axis.  (This  nonrobustness  is 
not  crucial.)  As  above,  if  only  one  Og  value  was  originally  stipulated,  then  we  are 
not  asked  for  it  again  here. 

alpha.  We  also  enter  a  value  for  a  (degrees)  if  it’s  asked  for  (as  for  bearing  and  9s)- 

delta_tReceiver  values.  Here,  if  asked,  we  enter  the  time  At  taken  to  fly  the  receiver 
flight  path,  in  seconds.  Bearing  measurements  will  be  taken  at  one  second  intervals. 

After  entering  these  values  the  relevant  graph  is  plotted.  We  can  also  plot  relative  down- 
and  cross-range  errors  and  CEPs  as  functions  of  two  variables,  using  PlotMatrixSurf  aces  .m, 
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